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ABSTRACT 


•^iA  new  two  dimensional  panel  method  has  been  developed. 

This  method  uses  a  new  approximating  element;  the  circular 
arc,  and  a  new  singularity  representation;  the  sine  series, 
and  all  integrations  are  performed  analytically  for  maximum 
computational  efficiency.  The  method  was  applied  to  a 
circular  cylinder  and  to  several  different  types  of 
airfoils,  and  a  number  of  characteristics  which  define  the 
method  were  varied  to  determine  their  effects  on  the  ^ 
solution . 

The  body  is  represented  by  a  series  of  circular,  arcs  in 

which  are  defined  by  sets  of  three  points  on  the  surface. 

The  singularity  distribution  is  modeled  by  a  power  series 
expansion  in  terms  of  the  sine  of  an  angular  variable  which 
is  related  to  the  arc  length  of  each  panel.  The  method  was 
applied  to  the  problem  of  flow  over  a  circular  cylinder,  and 
characteristics  which  define  the  method  were  varied. 
Results  indicated  that  accuracy  was  not  significantly 
affected  by  the  type  of  singularity,  while  dramatic 
reductions  in  velocity  errors  were  achieved  by  increasing 
the  number  of  terms  in  the  singularity  series.  Further, 
increasing  the  number  of  panels  also  increased  the  accuracy 
of  the  solution,  the  effect  of  singularity  continuity  was 
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more  apparent  m  the  smoothness  of  the  resulting  velocity 
distributions  than  in  the  accuracy  of  the  solutions,  the 
method  was  not  critically  sensitive  to  control  point 
location,  and  the  method  was  found  to  be  computationally 
efficient  as  the  number  of  terms  in  the  series  was 
increased. 

The  method  was  then  applied  to  a  Joukowski  airfoil,  a 
NACA  0024  airfoil,  a  thin  symmetric  airfoil,  and  to  both  a 
symmetric  and  a  cambered  Karman-Tref ftz  airfoil. 

v 

- — • — -^Major  conclusions  from  this  study  were  that  the  method 


produced  very  accurate  solutions  over  the  major  part  of  the 
airfoil,  error  reduction  occured  as  both  the  number  of 
panels  and  the  number  of  terms  in  the  series  were  increased, 
the  effect  of  point  source  location  was  large  but  was  local 
and  could  be  controlled,  the  method  was  generally 
insensitive  to  minor  variations  in  panelling,  and  the 
accuracy  of  the  solution  increased  as  panel  curvature  was 
increased  from  relatively  flat  to  circular. 
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COMPUTATION  OF  INCOMPRESSIBLE  POTENTIAL 


FLOW  OVER  AN  AIRFOIL  USING  A  HIGH  ORDER 
AERODYNAMIC  PANEL  METHOD  BASED  ON 
CIRCULAR  ARC  PANELS 

I .  Introduction 

The  central  problem  in  aerodynamics  is  to  predict  the 
pressures,  forces,  and  moments  exerted  on  a  body  immersed  in 
a  flowing  fluid.  One  would  like  to  be  able  to  solve  the 
full  Navier-Stokes  equations  for  any  configuration  at  one's 
desk,  but  this  is  not  possible  today.  Fortunately  the  real 
needs  of  the  engineering,  research,  and  development 
community  both  in  the  Air  Force  and  in  industry  allow  this 
problem  to  be  approached  from  several  different  levels.  At 
one  level  is  the  engineer  who  requires  the  details  of  a  full 
viscous  solution,  and  is  willing  to  spend  the  time  and 
computer  resources  required,  and  to  accept  the  limitations 
in  geometric  complexity  which  in  some  cases  are  necessary, 
in  order  to  obtain  solutions  of  this  nature.  At  the  other 
extreme  is  the  engineer  involved  in  perhaps  a  preliminary 
design  application.  His  requirement  is  for  very  rapid 
solutions  for  general  configurations  which  can  be  used  to 
develop  airfoil  or  aircraft  performance  characteristics.  He 


might  also  require  the  capability  for  rapid  development  of 
parametric  studies  to  assess  effects  of  small  changes  in 
geometry  or  flight  conditions  on  the  flow  over  an  airfoil,  a 
wing,  or  a  full  aircraft  configuration. 

While  it  is  true  that  much  progress  has  been  made  over 
the  last  several  years  in  the  development  of  both  Navier- 
Stokes  solutions  and  non-linear  potential  solutions,  these 
areas  cannot  as  yet  satisfy  the  engineering  requirements 
described  above.  For  this  reason  much  interest  and 
attention  has  been  (and  continues  to  be)  focused  on  the 
development  and  improvement  of  linear  potential  flow 
solutions  in  general,  and  in  the  panel  method  approach  to 
obtaining  such  solutions  in  particular.  The  features  of  the 
panel  method  approach  which  make  it  particularly  attractive 
are  its  computational  efficiency,  and  its  ability  to 
accommodate  accurate  geometric  modeling.  In  addition,  it 
has  been  found  that  the  linear  potential  flow  model  provides 
sufficient  accuracy  for  many  engineering  applications,  and 
indeed  the  panel  method  approach  is  used  on  a  daily  basis  by 
industry  and  government  workers  to  solve  a  wide  variety  of 
two  and  three  dimensional  aerodynamic  problems. 

Given  the  unquestioned  value  and  utility  of  the  panel 
method  approach  to  solving  the  linear  potential  flow 
problem,  the  general  goal  of  this  dissertation,  which  will 
be  discussed  further  in  this  chapter,  is  to  develop  and 
investigate  a  particular  panel  method  approach  in  order  to 


add  to  the  level  of  understanding  of  such  methods.  The 
purpose  of  this  chapter,  then,  is  to  review  the  importance 
of  this  theory  and  discuss  the  assumptions  which  led  to  it; 
to  formulate  the  mathematical  statement  of  potential  flow 
over  a  body;  to  review  and  note  limitations  of  current 
methods  for  solving  this  problems;  and  to  present  the 
objective  of  this  dissertation. 

Linearized  Potential  Flow 

The  ability  of  the  linearized  potential  flow  equation 
to  accurately  model  flow  fields  about  realistic  flight 
vehicle  configurations  over  a  wide  range  of  realistic  flight 
vehicle  configurations  over  a  wide  range  of  flight 
conditions  is  well  known  to  aerodynamic ists.  These  results 
are  used  in  two  ways.  First,  they  can  predict  lift, 
moments,  and  induced  drag  for  complex  vehicles,  and  second, 
they  can  be  used  as  input  to  boundary  layer  calculations 
which  will  predict  friction  drag  and  separation.  In  fact, 
the  accuracy  of  boundary  layer  calculations  is  generally 
dependent  on  the  accuracy  of  the  input  potential  flow 
solution. 

The  basic  assumptions  leading  to  potential  flow  are 
that  the  fluid  is  inviscid,  non-heat  conducting,  isentropic, 
and  irrotational .  The  success  of  the  theory  lies  in  the 
fact  that  for  the  flow  of  air  over  a  body  the  effects  of 
viscosity  and  heat  transfer  are  confined  to  the  boundary 
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applicabl  e,  and  in  fact  form  the  basis  of  the  solution 
method  to  be  used  in  this  study. 

A  crucial  step  in  the  solution  of  Laplace’s  equation  is 
the  application  of  boundary  conditions.  In  small 
disturbance  theory  the  boundary  conditions  are  often  applied 
on  a  plane  rather  than  on  the  actual  surface  of  the 
configuration.  Although  this  linearization  of  the  boundary 
conditions  is  justified  if  the  body  is  thin  and  if  the  angle 
of  attack  is  small,  it  produces  a  non-physical  singularity 
at  the  airfoil  leading  edge.  This  singularity  can  be 
removed,  however,  by  applying  the  boundary  conditions  on  the 
actual  surface  of  the  configuration,  even  though  the 
governing  equation  was  derived  using  small  perturbation 
assumptions.  This  ability  to  apply  boundary  conditions 
exactly  and  in  a  convenient  way  is  an  important  feature  of 
the  panel  method  approach  to  solving  the  potential  flow 
problem. 

Statement  of  the  Problem 

The  mathematical  problem  of  potential  flow  about  the 
exterior  of  a  body  may  be  formulated  in  the  following 
manner.  Consider  a  closed  surface,  S,  (Figure  1)  immersed 
in'  a  flow  with  free  stream  velocity  ^  .  Let 

0O 

V  *  Vro  +  v  (3) 

where  f  is  the  total  velocity  and  v  is  a  perturbation 
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velocity.  Now  v  is  irrotational ;  therefore 

v  =  V4>  (4) 

where  <}>  is  the  perturbation  potential.  Assuming  .'1,-0 
for  convenience,  the  governing  equation  for  the  flow  becomes 

?2<J>=  0  (5) 

in  the  region  exterior  to  S  .  The  boundary  conditions  for 
this  problem  are  that  the  surface  is  impermeable,  and  the 
perturbation  velocities  are  zero  at  infinity.  That  is, 

|  V0 1  ■+■  0  at  infinity  (6a) 

V<J>»n  +  V'oo*n  =  OonS  (6b) 

where  n  is  the  outward  surface  normal  vector  on  S 
Once  a  solution  for  <P  is  determined,  the  pressure  on  S 
is  found  from  Bernoulli's  equation  as 


Many  approaches  have  been  used  to  solve  the  problem 
posed  in  equations  five  and  six,  including  conformal 
mapping,  finite  difference,  finite  element,  and  singularity 
methods.  Conformal  mapping  methods  (Refs  3,4)  have  been 
used  to  obtain  accurate  solutions  in  two  dimensional  cases, 


but  they  can  not  be  extended  to  three  dimensions.  The 
problem  with  using  a  conformal  mapping  approach  is  the 
generation  of  mappings  for  arbitrary  shapes. 

There  are  many  finite  difference  methods  (Refs  5, 6, 7, 8) 
which  solve  the  exterior  potential  flow  problem,  but  they 
are  usually  applied  to  the  linear  one  represented  by  Eq  5. 
These  methods  use  transformations  to  map  the  physical  space 
into  a  computational  space  in  which  boundary  conditions  can 
be  applied  with  less  difficulty.  Results  for  general  two 
dimensional  shapes  have  been  obtained,  but  the  methods  have 
only  limited  ability  to  handle  complex  three  dimensional 
geometry.  Disadvantages  of  finite  difference  methods  as 
applied  to  either  the  linear  or  non-linear  problem  are  that 
the  solution  must  be  found  throughout  the  entire  flow  field 
and  that  computer  time  and  storage  requirements  are  large 
(even  in  two  dimensions). 

A  newer  approach  is  the  finite  element  technique  (Refs 
9,10,11).  Developed  initially  as  a  structural  analysis 
tool,  there  has  been  considerable  application  of  the  method 
to  fluid  dynamics  problems.  It  has  had,  however,  relatively 
little  application  to  problems  involving  complex  geometric 
configurations,  and  shares  with  finite  difference  methods 
both  the  disadvantage  of  requiring  a  solution  throughout  the 
flow  field,  and  the  advantage  of  being  applicable  to  the 
non-linear  formulation  of  the  problem. 

Finally,  singularity  methods  have  been  used  for  many 
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a  plane  surface  (for  example,  the  camber  line  of  a  wing  or 
airfoil)  with  no  thickness  to  model  the  desired 
configuration. 

Compared  to  a  discrete  singularity  method,  a  panel 
method  (or  actually  any  distributed  singularity  method)  has 
advantages  which  are  related  to  the  order  of  the 
singularity.  One  cannot  compute  flow  quantities  at  a  point 
singularity  because  they  are  mathematically  undefined  there. 
One  can,  however,  make  such  computations  on  a  panel 
containing  a  distributed  singularity,  except  at  the 
endpoints  of  the  panel  where  the  flow  quantities  are  again 
singular.  However,  this  singularity  is  of  a  lower  order  than 
the  point  singularity.  This  means  that  for  a  given  level  of 
accuracy  one  can  perform  computations  closer  to  the 
endpoints  of  the  distributed  singularity  panel  than  one  can 
to  the  point  singularity.  This  is  important  because  most 
configurations  of  interest  contain  regions  where  one  part  of 
the  surface  is  near  another  part  of  the  surface.  An  example 
would  be  the  upper  and  lower  surfaces  at  the  trailing  edge 
of  a  wing  or  airfoil. 

Compared  to  finite  difference  (FD)  and  finite  element 
(FE)  methods,  panel  methods  may  require  less  computer 
resources  for  a  given  configuration  and  level  of  accuracy. 
The  reason  for  this  is  that  to  solve  a  three  dimensional 
problem,  the  FD  or  FE  method  must  solve  a  partial 
differential  equation  in  three  independent  variables. 
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This  requires  generation  of  a  mesh  of  grid  points  which 
fills  the  volume  of  the  flow  field.  To  solve  the  same 
problem  ,  the  panel  method  must  solve  a  two  dimensional 
integral  equation.  This  requires  generation  of  a  mesh  of 
grid  points  only  on  the  surface  of  the  configuration.  In 
effect,  the  dimension  of  the  problem  that  must  be  solved  has 
been  reduced  by  one.  This  reduction  also  occurs  if  the 
original  problem  is  that  once  the  panel  method  solution  has 
been  obtained,  the  flow  quantities  can  be  determined  at  any 
other  point  in  the  flow  field  by  a  simple  matrix 
multiplication.  This  is  in  contrast  to  the  FD  or  FE  methods 
for  which  the  solution  in  the  flowfield  is  computed 
certainty  only  at  the  grid  points  used  to  obtain  it.  One 
could,  of  course,  interpolate  these  values,  but  to  obtain 
velocities,  for  example,  from  a  solution  for  potential,  one 
would  have  to  use  a  numerical  differentiation  scheme  of  some 
sort  which  would  introduce  additional  inaccuracy  into  the 
result . 

Another  advantage  of  the  panel  method  compared  to  a  FD 
method  is  that  the  panel  method  Can  often  model  complicated 
geometries  more  easily.  The  reason  for  this  is  that  FD 
methods  often  require  a  coordinate  system  which  is  fitted  to 
the  configuration  surface  in  order  to  simplify  application 
of  boundary  conditions.  Generation  of  this  coordinate 


system  can  in  itself  require  the  numerical  solution  of  a  set 


of  partial  differential  equations.  The  panel 
method, however,  requires  only  the  surface  geometry  as  input. 

Although  the  advantages  of  panel  methods  as  described 
above  are  important,  it  must  be  remembered  that  the  method 
gives  a  linearized  potential  flow  solution  to  a  given 
problem.  Both  FD  and  FE  methods  are  applicable  to  the  non¬ 
linear  problem  (potential  and  non-potential)  as  well.  There 
are  many  situations  in  which  a  non-linear  solution  about  a 
simplified  configuration  is  more  useful  than  a  linearized 
potential  flow  solution  about  a  more  detailed  and  exact 
geometric  representation  of  the  configuration.  Conversely, 
it  is  also  true  that  there  are  a  great  many  applications  for 
which  the  linearized  potential  flow  solution  is  satisfactory 
and  in  these  cases  its  characteristics  of  geometric 
complexity  and  computational  efficiency  are  highly 
desireable . 

Literature  Review 

The  basic  theory  behind  the  use  of  a  panel  method  to 
solve  the  potential  flow  problem  (a  review  of  which  is 
given  by  Hess,  Ref.  12)  was  developed  from  the  principles  of 
potential  theory  (Ref.  13,14).  The  practical  application  of 
the  method  was  not  feasible,  however,  until  the  digital 
computer  became  available.  Since  the  early  1960 's  work  in 
this  area  has  increased  greatly  from  initial  efforts  at 
computing  axisymmetric  and  non-lifting  three  dimensional 


flow,  through  higher  order  lifting  two  dimensional  flow,  and 
lifting  three  dimensional  flow,  to  the  present  day  where 
complex  configurations  are  being  calculated  in  supersonic 


flow.  This  review  will  cover  the  major  two  dimensional 
panel  methods  available  today,  followed  by  a  discussion  of 
representative  three  dimensional  work,  and  will  conclude 
with  a  summary  of  limitations  in  the  methods  available  at 
the  present  time. 

Two  Dimensional  Methods.  A  large  number  of  two 
dimensional  methods  have  been  developed  over  the  past 
several  years.  Since  many  of  them  are  similar  in  concept,  a 
representative  sample  illustrative  of  different  approaches 
has  been  selected  for  discussion.  Table  I  presents  some 
general  characteristics  and  unique  features  of  these 
methods. 

Hess's  low  order  method  (Ref  18)  used  constant  source 
and  vorticity  distributions  on  flat  panels.  His  higher 
order  method  (Ref  15,16)  models  the  surface  as  a  series 
expansion,  truncated  such  that  the  representation  is 
parabolic,  while  the  integrand  in  the  velocity  influence 
integral  is  expanded  in  a  series  that  assumes  the  surface  is 
nearly  flat.  The  method  also  used  source  and  vorticity 
distributions  in  which  the  vorticity  is  assumed  to  vary 
parabolically  in  arc  length  from  the  trailing  edge  of  the 
airfoil  through  the  leading  edge  and  back  to  the  trailing 
edge  where  it  is  zero.  The  higher  order  method  shows 
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Vorticity  Representa¬ 
tion. 


increased  accuracy  over  the  lower  order  method,  particularly 
for  internal  flows. 

Henshaw  (Refs  19-28)  has  developed  a  variant  of  Hess's 
higher  order  method  which  used  a  quartic  polynomial 
representation  of  the  surface,  and  which  expands  the 
velocity  influence  integrand  about  a  circular  arc  rather 
than  about  a  flat  panel.  He  reports  an  improved  accuracy 
with  this  formulation  but  his  results  are  difficult  to 
interpret.  Henshaw  has  also  formulated  an  approach  using 
vorticity  only,  with  an  error  parameter  which  allows 
specification  of  circulation.  This  parameter  is  added  to 
the  left  hand  side  of  the  boundary  condition  equations  with 
a  coefficient  which  is  specified  according  to  certain 
criteria. 

Bristow  (Refs  29-31),  using  Hess's  basic  as  well  as  his 
higher  order  method  has  formulated  two  interesting 
approaches  to  the  problem.  Using  the  basic  method,  he  has 
incorporated  a  singularity  strength  minimization  procedure 
which  reduces  source  strength  gradients,  and  thus  errors  in 
tangential  velocities.  Using  the  higher  order  method,  his 
formulation  allows  a  priori  determination  of  the  source 
strengths,  coupled  with  an  error  parameter  approach  to 
obtaining  the  vorticity  strengths.  The  second  method  also 
produces  singularity  strengths  with  mild  gradients  and  good 
accuracy,  but  at  lower  computing  cost  than  the  first.  Both 
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of  these  methods  have  a  design  capability  as  well. 
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The  methods  of  Raj  (Ref  32)  and  Keller  (Ref  34)  provide 
different  approaches.  Raj  used  a  piecewise  linear  vorticity 
distribution  on  a  surface  described  numerically,  and  all  the 
integrations  are  performed  numerically  in  the  physical 
airfoil  plane.  His  results  are  accurate,  but  the  method  is 
time  consuming.  Keller's  approach  is  to  generate  a 
transformation  which  maps  an  airfoil  into  a  near  circle.  He 
then  performs  all  integrations  numerically  in  the  circle 
plane.  This  is  advantageous  because  the  integrals  are 
easier  to  integrate  numerically  on  a  circular  or  nearly 
circular  surface  than  on  an  arbitrary  surface.  This  method 
is  not,  however,  extendable  to  three  dimensions. 

Three  Dimensional  Methods.  A  significant  amount  of 
work  has  been  done  in  the  area  of  three  dimensional  panel 
methods.  Characteristics  of  the  more  important  of  these  are 
shown  in  Table  II.  These  methods  will  be  discussed  further 
in  the  following. 

Hess's  method  (Refs  18,34,35)  was  the  first  surface 
paneling  method  applicable  to  arbitrary  geometries.  It  is 
an  incompressible  method  which  uses  flat  panels  to  model  a 
configuration.  Constant  sources  are  used  on  body  panels 
while  constant  vorticity  is  used  on  wing  panels.  The  wing 
panels  are  lumped  into  chordwise  strips  over  which  a 
parabolic  distribution  of  vorticity  is  placed,  so  that  only 
one  vorticity  unknown  is  associated  with  each  strip.  A 
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Kutta  condition  is  applied  at  the  trailing  edge  of  each 
strip  to  obtain  this  unknown,  while  the  source  unknowns  are 
found  by  applying  a  normal  velocity  boundary  condition  on 
each  panel.  Hess  has  also  developed  a  higher  order  three 
dimensional  method  (Ref  36)  which  is  basically  an  extension 
of  his  two  dimensional  work.  The  new  method  has  shown  an 
improved  accuracy  for  a  given  number  of  panels,  but  is  at 
present  a  non-lifting  method. 

Woodward's  basic  method  (Ref  37)  was  the  first  unified 
(subsonic  and  supersonic)  method  for  general  configurations. 
It  modeled  a  surface  with  flat  panels,  which  were  not 
necessarily  contiguous.  Linear  sources  and  constant 
vorticity  were  distributed  on  the  panels  to  account  for  lift 
effects  and  line  sources  and  doublets  provided  body 
thickness  effects.  A  normal  velocity  boundary  condition  was 
applied  at  a  control  point  whose  location  was  chosen  so  as 
to  provide  the  best  results.  This  method  was  successful, 
but  was  limited  in  the  degree  of  geometric  complexity  that 
it  could  model,  and  was  sensitive  to  control  point 
placement . 

In  1973  an  improved  version  of  the  method  (Ref  38)  was 
presented  which  retained  a  flat  panel  surface  representation 
but  used  linear  source  and  linear  vorticity  singularities, 
and  which  had  planar  and  non-planar  boundary  condition 
options.  Linear  sources  were  distributed  on  body  surface 


panels  for  both  options,  while  on  wing  panels  the  planar 
option  used  linear  source  and  vorticity  distributions,  and 


the  non-planar  option  used  only  linear  vorticity.  A  normal 
velocity  boundary  condition  was  applied  at  a  control  point 
located  at  a  panel  centroid.  This  method  allowed  more 
accurate  modeling  of  body  shapes  and  exhibited  reduced 
sensitivity  to  control  point  location.  There  was  a 
difficulty,  however,  in  using  the  non-planar  option  in 
supersonic  flow  because  the  panels  exhibited  discontinuities 
in  slope  and  position.  This  caused  disturbances  to 
propagate  downstream  inside  the  configuration  being  modeled 
(that  is,  in  the  non-physical  interior  flow)  in  such  a  way 
as  to  eventually  destroy  the  solution  on  the  exterior  of  the 
body. 

Woodward  has  recently  developed  a  solution  to  this 
problem  (Ref  39)  using  a  combined  source  and  vortex  called  a 
triplet  singularity.  This  singularity  controls  the  interior 
flow  by  cancelling  perturbation  velocities  there  without 
explicitly  applying  boundary  conditions  in  the  interior 
region.  This  approach  has  shown  good  results  when  applied 
to  bodies,  but  has  yet  to  be  applied  to  wings,  or  more 
general  configurations. 

Robert's  method  (Refs  40-42)  uses  surface  sources  and 
internal  doublet  sheets  to  compute  subsonic  flow  about 
general  configurations.  The  surface  is  mapped  to  a 


parametric  plane  where  it  is  represented  as  a  bicubic 


spline,  and  the  singularity  distributions  are  modeled  as 
bicubic  splines  on  this  surface.  This  approach  is  capable 
of  yielding  very  accurate  solutions,  but  all  the  mappings 
and  the  integrations  of  the  influence  coefficients  are  done 
numerically;  thus  an  extremely  large  amount  of  computer  time 
is  required  to  obtain  a  solution.  For  this  reason  the 
method  has  not  been  widely  used. 

Torino's  method  (Refs  43-45)  is  a  general  method  for 
unsteady  subsonic  or  supersonic  flow  over  arbitrary 
configurations.  It  uses  constant  source  and  doublet 
distributions  on  hyperboloidal  surface  panels,  with  an 
interior  potential  boundary  condition.  Preliminary  results 
using  this  method  seem  to  be  good,  but  it  has  not  been 
tested  extensively  to  date.  It  should  be  emphasized  that 
the  method  was  developed  to  solve  the  general  unsteady 
problem,  and  is  perhaps  the  most  advanced  in  this  area. 

Over  the  past  ten  years,  researchers  at  the  Boeing 
Company  have  developed  a  general  subsonic  and  supersonic 
method  applicable  to  arbitrary  configurations.  The  method 
has  evolved  from  a  low  order  subsonic  method  to  a  higher 
order  supersonic  method  known  as  the  PANAIR  (Paneling 
Aerodynamics)  system.  In  1967  Rubbert  et  al  (Ref  46) 
described  a  subsonic  method  using  flat  panels  with  a  surface 
distribution  of  constant  sources,  and  an  interior  doublet 
distribution.  The  method  produced  good  results,  but  the  use 
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of  the  internal  lifting  system,  coupled  with  the  use  of  flat 
panels  limited  the  degree  of  geometric  complexity  which 
could  be  easily  and  accurately  modeled. 

In  1972  Rubbert  and  Saaris  (Ref  47)  presented 
additional  results  using  the  same  basic  method,  but  with  the 
addition  of  internal  singularity  sheets  which  were  used  to 
maintain  an  internal  flow  which  (although  of  no  physical 
interest)  would  improve  the  external  flow  characteristics. 
This  method  was  sensitive  to  the  paneling  arrangement  since 
it  was  a  low  order  method. 

To  correct  some  of  these  problems,  Johnson  and  Rubbert 
(Ref  48)  developed  a  higher  order  subsonic  method.  Key 
features  were  the  use  of  linear  sources  and  quadratic 
doublets  distributed  on  curved  panels,  with  internal  doublet 
sheets  to  provide  lift  effects.  Since  these  panels  were 
developed  by  fitting,  in  a  least  squares  sense,  parabolic 
curves  through  the  actual  surface  points,  the  panels  were 
discontinuous.  Further,  the  panels  were  restricted  to  being 
only  slightly  curved  through  the  use  of  a  near  field 
expansion  for  calculating  the  influence  coefficients.  A  far 
field  expansion  was  also  used  to  increase  computational 
efficiency.  Results  were  obtained  for  a  randomly  paneled 
sphere  and  wing,  which  indicated  the  versatility  of  the 
method. 

The  method  was  then  extended  to  supersonic  flow  by 
Ehlers,  Johnson,  and  Rubbert  (Ref  49)  using  linear  sources 
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and  quadratic  doublets  on  slightly  curved  panels  with  a 
linearized  mass  flux  boundary  condition  instead  of  the  usual 
velocity  boundary  condition.  In  addition,  an  internal 
potential  boundary  condition  was  used  to  control 
disturbances  in  the  interior  flow  which  tended  to  amplify  as 
they  were  reflected  by  the  interior  surface  and  which  would 
perturb  or  destroy  the  exterior  solution.  This  approach 
also  allowed  a-priori  determination  of  source  strengths 
which  reduced  the  order  of  the  system  of  linear  equations 
which  had  to  be  solved.  Results  were  presented  for  a 
randomly  paneled  spindle,  an  inlet  with  nacelle,  and  several 
wings  which  showed  excellent  agreement  between  experiment  or 
theory,  and  the  computed  results.  A  problem  developed 
however,  because  both  the  paneling  and  the  doublet 
distribution  (which  was  found  using  a  least  squares 
approach)  were  slightly  discontinuous.  This  generated 
singularities  which  propagated  along  Mach  lines  with 
undiminished  strength,  and  which,  if  downstream  control 
points  were  too  close,  could  cause  the  influence  coefficient 
matrix  to  be  singular. 

The  solution  to  this  problem,  given  by  Ehlers  et  al 
(Refs  50-52)  was  to  replace  the  discontinuous  curved  panel 
concept  by  a  continuous  flat  panel  concept.  Previous  flat 
panel  methods  used  four  input  corner  points  to  define  a 
single  flat  panel  which  did  not  necessarily  pass  through  the 


input  points.  This  new  method  used  four  input  corner  points 
to  define  five  planar  subpanels  which  passed  through  the 
corner  points  and  which  were  continuous  with  all  neighboring 
panels.  The  result  modeled  a  surface  with  continuous  flat 
panels,  and  allowed  the  quadratic  doublet  distributions  on 
each  subpanel  to  be  exactly  continuous  along  all  edges. 
This  method  has  produced  good  results  to  date  (Ref  53),  and 
is  the  first  higher  order  supersonic  method  capable  of 
accurately  modeling  extremely  complex  geometric 
configurations. 

In  general,  low  order  three  dimensional  panel  methods 
are  fairly  complex,  require  large  numbers  of  panels  to 
achieve  a  reasonable  accuracy,  and  are  sensitive  to  panel 
and  control  point  placement.  The  higher  order  methods  have 
reduced  the  number  of  panels  required  to  achieve  the  same  or 
better  accuracy,  but  at  a  cost  of  increased  complexity  and 
computational  requirements.  Some  comparisons  between 
several  of  these  methods  are  given  by  Thomas  and  Miller  (Ref 
54),  and  Landrum  and  Miller  (Ref  55). 

Objective  of  Dissertation 

The  motivation  for  studying  linear  potential  flow 
methods  in  general  and  panel  methods  in  particular  stems 
from  the  proven  usefulness  of  these  methods  in  a  wide  range 
of  engineering  and  research  activity.  It  has  also  been 


shown  that  current  available  methods  are  in  general 


complicated  require  in  many  cases  large  numbers  of  panels  to 
achieve  a  given  level  of  accuracy,  and  require  a  large 
degree  of  expertise  on  the  part  of  the  user  in  order  to 
obtain  satisfactory  results.  A  partial  reason  for  these 
deficiencies  is  that  the  influence  of  a  number  of  the 
characteristics  that  define  the  panel  method  approach,  in 
both  two  and  three  dimensional  cases,  are  not  adequately 
understood.  These  characteristics  include  panel  curvature, 
singularity  distribution  continuity,  type  of  singularity, 
order  of  the  singularity  approximation,  the  type  of  boundary 
condition,  the  numerical  implementation  of  the  Kutta 
condition,  and  control  point  location.  The  question  of  the 
effects  of  these  characteristics  on  solution  accuracy  for  a 
given  method  has  not  been  fully  answered. 

The  objective  of  this  dissertation  is  to  answer  these 
questions  within  the  framework  of  a  two  dimensional 
incompressible  method  as  a  first  step  in  developing  a  fuller 
understanding  of  the  effects  of  these  characteristics.  The 
results  of  such  an  investigation  will  provide  guidance  to 
others  who  wish  to  develop  two  or  three  dimensional  panel 
methods  for  their  own  specific  applications.  To  accomplish 
this  objective,  a  new  two  dimensional  method,  based  on  the 
use  of  circular  arc  panels,  has  been  formulated,  and  has 
been  extensively  tested  in  applications  to  the  cases  of  flow 
over  a  circular  cylinder  and  flow  over  several  types  of 


airfoils. 


The  results  have  shown  that  accuracy  increased  as 
additional  terms  in  the  series  representing  the  singularity 
distribution  are  kept,  as  panel  curvature  is  varied  from 
flat  to  circular  and  as  continuity  of  the  vorticity 
distribution  is  enforced.  Additionally,  the  effect  of 
control  point  location  has  been  found  to  be  relatively 
small,  and  the  required  number  of  panels  for  a  given 
accuracy  has  been  found  to  be  less  than  that  required  by  the 
method  of  Raj.  The  present  method  has  also  been  compared  to 
Hess's  higher  order  method  as  formulated  by  Bristow  for  a 
thin  airfoil,  and  has  been  found  to  give  a  small  improvement 
in  computed  perturbation  velocities. 

The  following  chapters  of  this  dissertation  will 
discuss  these  points  in  detail.  The  next  chapter  will 
briefly  present  some  highlights  of  potential  theory,  and  the 
general  panel  method  approach.  Then  the  details  of  a  new 
paneling  method  based  on  circular  arc  panels  will  be 
presented,  followed  by  the  application  of  this  method  to  the 
circular  cylinder  problem,  and  then  to  the  airfoil  problem. 
Finally  the  conclusions  resulting  from  this  work  and  ideas 
concerning  possible  extensions  of  the  method  will  be 
presented. 
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II .  Two  Dimensional  Potential  Theory 


The  purpose  of  this  chapter  is  to  present  briefly  some 
basic  facts  about  potential  theory  and  surface  singularity 
distributions  which  will  have  direct  application  to  the  work 
that  follows.  It  will  be  seen  that  the  panel  method 
approach  to  solving  flow  problems  is  dependent  on  the 
results  of  potential  theory.  Harmonic  functions  will  first 
be  discussed  from  a  partial  differential  equation  viewpoint, 
followed  by  a  presentation  of  the  Green's  theorem 
representation  of  a  harmonic  function.  Some  characteristics 
of  surface  singularity  distributions  and  some  relevant 
properties  of  integral  equations  will  then  be  discussed. 
Finally  the  reduction  of  the  Green's  theorem  formulation  to 
an  integral  equation  will  be  considered  along  with  some 
unresolved  questions  which  arise  in  this  formulation. 

Harmonic  Functions 

Solutions  to  Laplace's  equation  are  called  harmonic 
functions.  Such  functions  have  properties  which  allow  the 
development  of  the  integral  equation  method  which  is  the 
basis  of  the  panel  method  approach  to  solving  Laplace's 
equation.  This  equation  has  also  been  studied  considerably 
from  a  partial  differential  equation  viewpoint  in  which  one 
determines  conditions  which  guarantee  the  existence  and 


uniqueness  of  solutions  to  a  given  equation  which  is  subject 
to  a  particular  specification  of  boundary  conditions. 

The  general  boundary  value  problem  can  be  stated  as 
follows:  find  a  function  <J>  which  satisfies  V 2 0  =  0  in  a 
region  R  and  where  either  4>  =  f(s)  or  =  f(s)  on  the 
boundary  of  R,  and  where  f(s)  is  a  known  function.  If 

is  specified,  this  problem  is  called  a  Dirichlet 
problem;  while  if  -yj.  is  specified,  it  is  called  a  Neumann 
problem.  The  existence  and  uniqueness  of  solutions  to  these 
problems  depends  on  whether  R  is  an  interior  or  exterior 
region.  Given  that  the  boundary  values  are  continuous  and 
that  the  boundary  is  sufficiently  regular,  Table  III  (Ref 
14)  summarizes  the  conditions  for  which  these  problems  have 
solutions. 

These  results  from  the  theory  of  partial  differential 
equations  will  be  used  to  verify  the  correctness  of  the 
integral  equation  formulation  of  the  problem  which  leads  to 
the  panel  method  solution  to  Laplace's  equation.  This 
formulation  is  dependent  on  the  property  of  harmonic 
functions  that  is  state!  in  Green's  theorem. 

Green's  Theorem  Formulation 

Using  Green's  theorem,  the  value  of  a  harmonic  function 
at  any  point  in  a  region  may  be  expressed  in  terms  of  its 
value  on  the  boundary  of  the  region.  This  form  may  then  be 
interpreted  as  a  singularity  distribution  on  the  boundary. 
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TABLE  III 


Boundary  Value  Problem  Solutions 


Boundary 

Value  Prob. 

Region 

Interior  Exterior 

Dirichlet 

Solution 

No  Solution 

Neumann 

Solution* 

Solution 

*  If  and  only  if  Jr  dl  =  0 

B 

Green's  theorem  may  be  applied  to  harmonic  functions  which 
are  single  valued  in  some  region.  The  problem  may  be 
formulated  in  two  dimensions  by  considering  (Figure  2)  two 
harmonic  functions,  $  and  <j>  j  ,  and  two  regions,  R  and 
Ri  which  are  divided  by  a  boundary  curve  B.  The  curve  b  is 
called  a  barrier  and  is  required  to  make  R  a  simply 
connected  region  which  then  insures  that  <{>  will  be  single 
valued  there.  If  P  is  a  point  in  R,  it  can  be  shown  (Refs 
13,14,56)  that 

*(p)  =  ir  /(log  r  -  *!n  (log  r))  dl 


!_  f 
/ 


(4>+  -  4>”)  (log  r)  db 


Region  R' 
Potential  <j>  1 


00 


1  / 


Curve  B 


Figure  2.  Green's  Theorem  Formulation. 


and  also  that 


4>i(P)  -  ^  /  (log  r  -jyp-  -  <J>  i  |^(log  r )  )dl  =  0 


Adding  Eqs  8  and  9  gives 


<p(P)  -  — =•  f  (log  r(||'  -  -  (0-  $i)i_(i0g  r)  )dl 

B  3n 


(10) 


h  f  (*+  ~  <j>")  Tn  (log  r) 


Equations  8  and  10  show  that  <P  is  not  uniquely  determined 
until  both  <f>  and  <}>i  are  specified  on  B.  This  means  that 
the  solutions  for  R  and  Ri  are  independent  in  the  sense  that 
the  solution  could  be  changed  in  one  region  without  changing 
the  solution  in  the  other  region. 

Velocity  Potential.  If  <p  and  4> x  are  assumed  to  be 
velocity  potential  functions  and  if  a  =  and 


M  =  <P  ~  <f>l 


then  Eq  10  becomes 


4>(P) 


4?  /  (« 


log  r  -  y  |tr(log  r))  dl 


(11) 


+  h  f  (*+  -  *“>4r(1°e  r> 


The  first  integral  in  this  equation  can  be  interpreted  as 


potential  due  to 


source  distribution  of 


strength 


and  a  doublet  distribution  of  strength  y  on  B 
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while  the  second  can  be  interpreted  as  the  potential  due  to 


a  doublet  distribution  of  strength  y  on  b.  As  noted 
above,  <j>  is  not  unique  unless  both  a  and  y  are 
determined.  One  procedure  would  be  to  specify  a-priori 
either  a  or  y  ,  and  then  apply  another  boundary 
condition  to  determine  the  remaining  unknown.  This  is 
equivalent  to  specifying  the  solution  in  Rx  and  then  solving 
for  the  solution  in  R. 

Suppose  a  is  specified  on  B  to  be  a  =  0  .  Then  Eq 
10  becomes 


4>(P) 


1_ 

2tr 


/  n  in  (1°e  r)dl  +  §7  I  A*(|n(log  r))db 

B  b 

(12) 


where 


Ab 


<p 


Since  <J>  is  a  velocity  potential  A<j>  is  the  circulation 
around  B,  and  is  constant.  Also,  since  the  location  of  b  is 
arbitrary,  the  normal  and  tangential  derivatives  of  <j>  are 
continuous  across  b  (Ref  13).  Therefore  A <p  is  constant  on 
b,  and  the  second  integral  in  Eq  12  represents  a  constant 
strength  doublet  sheet  extending  to  infinity  and  is 
equivalent  to  a  wake.  If,  however,  y=0  on  B,  then  Eq  10 


becomes 


(13) 


<J>(P)  =  2tt  J  a  log  r  dl 

B 

In  this  case  y  =  0  implies  that  <j>  =  ^ x  on  B  and  thus 
that  <}>+  =  <p~  because  $  x  is  single  valued  in  Rj  ; 
therefore  the  integral  over  b  is  zero.  This  also  means  that 
one  cannot  obtain  circulation,  or  lift,  using  a  source 
distribution  only. 

Stream  Function.  In  the  last  section  <P  was  assumed 
to  be  a  velocity  potential  function,  although  the  general 
formulation  is  not  dependent  on  this  interpretation. 
Since  <P  may  be  any  harmonic  function,  assume  that  it 
represents  a  stream  function  ^  .  If  there  are  no  sources 
inside  B  then  ip  will  be  constant  (single  valued)  on  B 
because  it  is  a  measure  of  the  mass  flux  across  the  curve  B. 
Thus  ip+  =  ip  ,  and  the  expression  equivalent  to  Eq  10  is 


lp(p)  ~  w  /  (log  r  (!£  -  r>>  dl 

B 


(14) 


If  it  is  assumed  that  >p  -  ipi  on  B,  and  that 

_  dip  dipi 

r  “  "5 - JrT  then 


ip(P) 


-1 

Otf 


f 


Y  log  r  dl 


(15) 


Thus  the  stream  function  can  be  represented  in  terms  of  a 
surface  distribution  of  vorticity  with  strength  y 
Although  Eq  15  does  not  contain  a  wake  term  as  does  Eq  12, 
it  does  not  necessarily  represent  a  zero  circulation  case. 

Doublet-Vort icity  Sheet  Equivalence 

Hess  (Ref  35)  has  shown  that  the  velocity  field  due  to 
a  surface  distribution  of  doublets  (whose  axes  are  normal  to 
the  surface)  is  equivalent  to  the  combined  fields  of  a 
distribution  of  vorticity  on  the  surface  where  y  -  n  x  vu  ^ 
and  a  line  vortex  on  the  bounding  curve  of  the  surface  whose 
strength  is  equal  to  the  strength  of  the  doublet 
distribution  on  the  curve.  For  the  two  dimensional  case 
(Figure  3)  a  constant  doublet  sheet  of  strength  h  from  A 
to  B  is  equivalent  to  two  point  vortices  at  A  and  B  of 
strength  u  .  This  means  that  a  velocity  field  represented 
by  a  distribution  of  sources  and/or  doublets  can  also  be 
represented  by  a  distribution  of  sources  and/or  vortices. 
In  the  case  of  a  constant  doublet  distribution  on  a  wake, 
the  equivalent  vorticity  distribution  consists  of  a  pair  of 
point  vortices,  one  at  the  start  of  the  wake,  and  one  at 
infinity.  On  the  wake  y  “constant  implies  ?y  =  0  ,  and 

thus  y  *  0  on  the  wake.  Therefore  Eqs  12  and  15  are 
consistent  and  in  fact  are  equivalent. 
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Figure  3.  Doublet /Vortex  Sheet  Equivalence 


Singularity  Behavior 

Using  Green's  theorem  it  is  clear  that  the  problem  of 
potential  flow  over  a  body  can  be  modeled  using  several 
types  of  singularity  distributions.  These  surface 
distributions  exhibit  certain  properties  which  affect  how 
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Source  Sheet 


Vortex  Sheet' 


Figure  4  -  Properties  of  Singularity  Sheets 


they  may  be  used  to  model  different  types  of  flow. 

Transfer  of  boundary  conditions.  Consider  the 
singularity  distributions  shown  in  Figure  4  where  subscripts 
e  and  i  stand  for  surface  exterior  and  interior, 
respectively.  For  the  source  sheet,  the  potential  and  the 
tangential  velocity  are  continuous  across  the  sheet,  while 
the  normal  velocity  is  discontinuous.  For  the  vortex  sheet 
the  opposite  is  true;  that  is,  the  potential  and  the 
tangential  velocity  are  discontinuous  across  the  sheet  while 
the  normal  velocity  is  continuous.  The  importance  of  these 
properties,  as  emphasized  by  Rubbert  (Ref  57)  is  that  they 
cause  certain  characteristics  to  be  transferred  across  the 
surface.  For  example,  consider  a  source  sheet  on  a  closed 


surface.  If  somehow  a  distribution  is  specified  that  gives 
a  particular  solution  in  the  exterior  region,  the  resultant 
exterior  tangential  velocity  distribution  is  carried  across 
the  sheet  and  becomes  the  boundary  condition  specification 
on  the  interior  region.  In  the  case  of  a  vortex 
distribution,  the  normal  velocity  is  transferred  across  the 
sheet  so  that  the  interior  problem  becomes  effectively  a 
Neumann  problem.  But  recall  that  the  condition  ensuring  a 
solution  to  this  problem  is  that  the  net  normal  velocity,  or 
flow,  be  zero,  which  is  simply  a  statement  that  an 
incompressible  fluid  cannot  be  pumped  into  a  closed  region. 
One  procedure  for  alleviating  the  problem  of  a  non  zero  net 
normal  flow  would  be  to  place  a  sink  inside  the  surface  to 
remove  any  excess  fluid. 

Singularity  Behavior  at  Corners.  The  above  properties 
of  singularity  sheets  apply  to  surfaces  which  are  smooth  to 
some  order.  However,  many  bodies  of  interest  have  slope 
discontinuities  at  one  or  more  points,  such  as  an  airfoil 
with  a  sharp  trailing  edge.  Craggs  and  Mangier  (Ref  58) 
have  studied  the  behavior  of  source  distributions  at  corner 
points.  They  find  that  the  source  distribution  behaves  as  a 
power  of  distance  to  the  point  with  the  value  of  the  power 
depending  on  whether  the  flow  is  symmetric  about  the  corner 
and  whether  the  corner  is  concave  or  convex  to  the  flow. 
For  the  case  shown  in  Figure  5,  which  is  symmetrical  flow 


Source  Sheet,  Strength  a 
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Figure  5.  Behavior  of  a  Source  Sheet  at  a  Corner 


about  a  convex  corner, 
that  c  -*■  0  as  s  -*■  0 


the  power  is  positive,  so 
Thus  in  modeling  an  airfoil  with 


a  source  distribution  this  behavior  should  be  considered. 

Integral  Equations  -  Existence  and  Uniqueness 

In  this  section  some  results  from  the  theory  of 
integral  equations  (Refs  59,  60,  61),  which  will  be  applied 
in  the  following  sections,  will  be  discussed.  Equation  16 
is  the  general  form  of  a  Fredholm  integral  equation  of  the 
second  kind  where  K(x,y)  is  a  given  kernel  function, 

b 

$ ( x )  -  xj K(x,y)<j>(y)  dy  =  f(x)  (16) 

a 

f(x)  is  a  given  function,  X  is  a  parameter,  and  <|>(x)  is 
the  unknown.  Several  results  can  be  stated  about  this 
equation. 

1.  Either  Eq  16  has  a  nontrivial  solution,  or  the 
associated  homogeneous  equation 

b 

w(x)  - \j  K(x,y)w(y)dy  =  0  (17) 

a 

has  a  nontrivial  solution.  The  values  of  X  for  which  Eq 
17  has  nontrivial  solutions  are  called  eigenvalues,  and  the 
solutions  w(x)  are  called  eigenfunctions. 

2.  If  X  is  an  eigenvalue  of  Eq  16  then  this  equation 


is  an  inconsistent  equation  (i.e.  has  no  solution)  unless 


b 

J u(x)  f (X)  dx  =  0 
a 


(18) 


where 


b 

u( x )  -A J' K(y,x)u(y)  dy  =  0 
a 


(19) 


3.  If  Eq  18  holds,  then  there  are  an  infinite  number 
of  solutions  to  Eq  16  of  the  form 


<f>(x)  =  <|> (x)  +  c  w  (x)  (2°) 

p  m  m  ' 

m 

where  <p  is  a  particular  soluti  n,  the  c  are  arbitrary 
P  m 

constants,  and  the  summation  extends  over  the  set  of 
linearly  independent  eigenfunctions,  .  Another 

important  property  of  a  Fredholm  equation  of  the  second  kind 
is  that  it  is  equivalent  to  a  system  of  linear  algebraic 
equations. 


While  Fredholm  integral  equations  of  the  second  kind 
have  some  very  nice  properties,  Fredholm  equations  of  the 
first  kind,  of  which  Eq  21  is  the  general  form,  do  not. 


b 

X  / K(x,y)$(y)dy  =  f(x)  (21) 

a 

f(x)  *  known  function 

It  can  be  shown  that  equations  of  this  type  do  not  always 


have  solutions,  and  the  solutions  of  solvable  cases  are 
often  not  unique.  This  question  can  be  related  to  the 
properties  of  the  Dirichlet  and  Neumann  problems  which 
required,  essentially,  that  there  be  no  net  flux  into  a 
closed  region.  Equations  16  and  21  will  be  used,  with  some 
modification,  to  solve  several  problems  in  the  succeeding 
sections . 

Reduction  of  the  Singularity  Distribution  Formulation  to 
an  Integral  Equation 

Consider  the  problem  of  two  dimensional  incompressible 
flow  about  a  body,  B,  immersed  in  a  free  stream,  ^  ,  as 

GO 

shown  in  Figure  6,  where  P  is  a  field  point,  q  is  any  point 

on  B,  r (P, q)  is  the  distance  between  P  and  q,  and  n  is 

Q 

the  outward  normal  to  B  at  q. 

Surface  Source  Distribution.  The  velocity  potential 
function  may  be  represented  as  a  source  distribution  on  B  by 


<kp>  =  Sr  /c(q)Ks(p,q)dl  (22) 

B 

where  K  (P,q)  is  the  source  kernel  function  in  two 
dimensions 

Ks(P,q)  =  log(r(P, q) ) 


It  can  be  shown  (Ref  14)  that  if  Eq  22  is  differentiated  and 


if  tiie  field  point  P  goes  to  a  surface  point  p  one  obtains 


■>  -  ^  -  h  /“«D^-(Ks(p,q))dl 


(23) 


where  3/3np  means  the  derivative  in  the  normal  direction 
to  B  at  p,  and  the  subscript  e  means  that  P  goes  to  p  in  the 
region  exterior  to  B.  Similarly, 


#  >.  *  ::2£Ei  -  h  /a<q>  dr(Ks(p-'1)tU 


(24) 


P  i 


total  velocity 


field 


is  ^  35  ^,+v  where  v  *  V<J>  ,  then  the  standard  boundary 
condition  of  zero  external  normal  velocity  can  be  written 


V  (p)  -  VM  (p)  +  v  (p)  *  0 
e  n  e 


(25) 


v  (p)  =  -Vw  (p) 
e  n 


(26) 


V(p)  *  In"  >  *  ^  -  57/0<q)T5-(KsCp-tl>)dl 


(27) 


therefore 
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o(p) 


1 

IT 


(28) 


/  ”<«)  4- 


B 


j|-(Ks(p,q))dl  -  -2V„  (p) 
P  n 


This  is  a  Fredholm  Equation  of  the  Second  Kind  for  the 
source  strength  a(p)  on  B.  Before  the  results  presented 
earlier  for  this  type  of  equation  can  be  applied,  the  kernel 


K(p,q)  =  (Ks(p,q)) 


must  be  considered.  A  cursory  examination  indicates 
that  K(p,q)  is  singular  at  q=p  ,  but  it  can  be  shown 
that  in  the  two  dimensional  case  the  singularity  is 
removable  if  the  curve  B  is  sufficiently  smooth. 

Conditions  for  Solvability.  What  constitutes 

sufficient  smoothness  is  not  completely  clear,  Tricomi  (Ref 
59)  and  Sternberg  and  Smith  (Ref  62)  specify  that  B  have 
continuous  curvature.  Mikhlin  (Ref  63)  and  Pogorzelski  (Ref 
64)  require  that  the  surface  satisfy  the  following 
conditions  (which  are  called  Liapunov  conditions): 

1.  The  surface  has  a  definite  normal  at  each  point. 

2.  There  exists  a  number  c>0  such  that  a  sphere  of 
radius  c  centered  at  a  point  on  the  surface  cuts  a  portion  of 
the  surface  such  that  every  line  parallel  to  the  normal  at 
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the  point  cuts  the  portion  only  once. 

3.  The  angle  between  the  normals  to  any  two  points  on 
the  surface  satisfies  the  following: 


8 |  £  arE 


where  8  is  the  angle  between  the  normals,  r  is  the 
distance  between  the  points,  and  a  and  e  are  positive 
constants.  This  question  will  arise  again  when  these 
theories  are  applied  to  shapes  of  actual  aerodynamic 
interest,  the  majority  of  which  have  at  least  one  point  of 
slope  discontinuity.  The  question  of  applicability  of  the 
theories  to  such  surfaces  has  not  been  satisfactorily 
resolved  to  the  author's  knowledge.  It  might  be  reasoned 
that  the  actual  viscosity  in  the  boundary  layer  will 
effectively  round  off  any  corners,  and  this  may  be  the 
answer.  Also,  although  all  methods  exhibit  decreased 
accuracy  in  the  trailing  edge  region,  the  quality  of 
solutions  over  the  remainder  of  the  airfoil  does  not  seem  to 
be  adversely  affected. 

If  it  is  assumed  that  B  has  the  requisite  smoothness, 
and  noting  that  A  =  j  is  not  an  eigenvalue  of  Eq  28,  it  is 
known  that  a  unique  solution  exists  for  any  given  free 
stream  flow.  The  problem  just  posed  is  equivalent  to  the 
Neumann  exterior  problem  seen  earlier  and  does  indeed  have  a 
unique  solution.  Recalling  the  discussion  of  the  properties 
of  source  sheets,  since  the  potential  is  continuous  across 


the  sheet,  a  Dirichlet  boundary  value  problem  is  effectively 
imposed  on  the  interior  of  B. 


Surface  Vorticity  Distribution.  Now  consider  the  same 
problem  assuming  a  vortex  distribution  on  B.  The  stream 
function  for  the  singularity  distribution  is 

<J>(P)  =  ^  f  Y(q)  log  ( r(P ,  q) )  dl  (29) 

B 

where  Y(q)  is  the  vortex  strength.  But  the  boundary 
condition  will  be  applied  in  a  way  first  suggested  by 
Martenson  (Ref  65).  Consider  the  total  stream  function  of 
the  flow, 


V  s  +  ip  (30) 

Now  on  a  streamline,  such  as  a  body  surface,  \p  =  constant, 
or 


-  n 
Tt  "  0 


(31) 


where  t  is  the  surface  tangent  direction.  Equation  31  is 
actually  a  statement  of  zero  external  normal  velocity.  To 
see  this,  consider  the  following  boundary  condition: 


which  is  a  statement  that  the  total  internal  tangential 


velocity  on  B  is  zero.  Now  Green's  theorem  for  ¥  states 
that 


//(,V2,+Vcc.V,)  ds  =  y* ¥  ||-  dl 


(33) 


where  S  is  interior  to  B.  Also, 


v2y  =0  in  S 


7H7?  =  V' 


where  V  is  the  total  velocity  inside  B.  Thus 


//v’ds  -/ 


*  Hr  dl 


(34) 


Now  the  boundary  condition  is 


3¥  _ 


=  0  on  B 


(35) 


so  that 


If 


V2ds  -  0 


(36) 


but  this  means  that  V  =  0  inside  B.  Further,  this  implies 
that  ^  ■  constant  inside  B,  and  thus 


4 


0  on  B 


(37) 


Now  is  the  interior  normal  velocity  on  B,  and  since 

1 

the  tangential  derivative  (normal  velocity)  is  continuous 
across  a  vortex  sheet 


on  B 


(38) 


This  is  just  a  statement  of  zero  exterior  normal  velocity, 
as  it  was  desired  to  show.  Now  computing  explicitly, 


f  Y(q)  r—  (log  r)  dl 
B  P 


and 


(39) 


liT  =  =1¥1  -  h  Y(q)  IT-  (1°e  r>  dl 

P  i  B  p 


(40) 


so  that  subtracting  these  gives 


*->  -i-)  -t(p) 


3n 


(41) 


Now  consider  the  total  external  tangential  velocity  given  by 
Eq  42  which  is  actually  the  quantity  of  interest. 
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Using  Eqs  32  and  41,  this  can  be  written 


1 

3 

S  9' 


y(p) 


(43) 


This  states  the  important  property  that  the  total  external 
tangential  velocity  is  equal  to  the  local  vorticity 
strength. 

Now  consider  the  integral  equation  results  as  they 
apply  to  this  formulation.  Applying  Eqs  32  and  40  one 
obtains 


ii) 


zl(,p),  _  I_  f 

2  2tt  J 


B 


—  9  ^ 

Y(q)  ^  (log  r)  dl  =  — 5~ 


y(p) 


+  I  /  v(q)  3K- 


(44) 


3T 

3np  (1°S  r>  dl  *  2  on 


This  is  again  a  Fredholm  equation  of  the  second  kind,  but 
now  the  parameter  A  =  ^  is  in  fact  an  eigenvalue.  From 
the  earlier  discussion  of  integral  equations  it  is  known 
that  when  a  solution  does  exist  ,  it  is  not  unique.  This 
non-uniqueness  will  be  removed  by  the  application  of  a  Kutta 
condition.  In  later  parts  of  this  work  a  vortex 


distribution  to  which  is  applied  the  standard  zero  exterior 


The  problem  of  two  dimensional  incompressible  potential 
flow  over  a  body  can  be  formulated,  using  the  concepts 
discussed  above,  in  terms  of  a  singularity  distribution  on 
the  surface  of  the  body.  The  different  singularity  types 
have  different  characteristics  which  determine  whether  they 
will  be  effective  in  a  particular  application.  Once  the 
singularity  has  been  chosen  and  the  problem  has  been  reduced 
to  the  appropriate  integral  equation,  additional  numerical 
approximations  must  be  introduced  in  order  to  obtain 
solutions  for  arbitrary  geometries.  The  details  of  these 
approximations  form  the  basis  of  the  panel  method  approach 
to  solving  this  problem. 

In  the  next  chapter  a  new  panel  method  will  be 
presented.  The  method  is  based  on  the  use  of  circular  arc 
panels  with  higher  order  singularity  distributions. 
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III.  Panel  Method  Approach 


The  purpose  of  this  chapter  is  to  formulate  a  new 
method  of  obtaining  an  approximate  solution  to  the  integral 
equations  developed  in  the  preceeding  chapter  using  the 
panel  method  approach.  This  method  is  based  on  the  concept 
of  approximating  the  surface  of  a  two  dimensional  body  by  a 
series  of  circular  arcs  on  which  higher  order  source  and 
vorticity  distributions  are  placed. 

Any  panel  method  consists  of  certain  assumptions  and 
approximations  concerning  the  basic  elements  of  the  integral 
equation.  These  elements  include  the  approximate 
representation  of  the  surface  over  which  the  integral  is 
taken,  the  approximate  representation  of  the  singularity 
distribution  which  is  assumed  on  the  surface,  the  type  of 
boundary  conditions  which  are  applied,  and  the  procedure  by 
which  the  Kutta  condition  is  satisfied. 

The  next  section  will  discuss  different  ways  of 
representing  the  surface,  and  will  give  the  rationale  for 
the  choice  of  the  circular  arc  element,  as  well  as  details 
of  the  numerical  implementation.  This  will  be  followed  by  a 
discussion  of  the  types  of  singularities  available  and  the 
procedure  by  which  the  singularity  strength  is  approximated 
on  a  panel.  A  discussion  of  the  boundary  conditions  that 


will  be  applied  will  follow.  This  will  include  the 
reduction  of  the  velocity  boundary  conditions  to  a  matrix 
equation,  and  a  discussion  of  the  form  of  the  Kutta 
condition  which  will  be  used.  Finally,  the  reduction  of  the 
method  to  a  system  of  linear  algebraic  equations  and  the 
procedure  by  which  the  system  is  solved  will  be  presented. 

Surface  Representation 

The  integral  equation  to  be  solved  contains  an  integral 
over  a  surface  for  which  an  analytic  description  will  not 
usually  be  available,  and  even  if  it  were  available,  the 
evaluation  in  closed  form  of  the  resulting  integral  would 
generally  be  impossible.  Therefore,  a  suitable 
approximation  to  the  surface  must  be  found  which  will 
capture  those  features  of  the  surface  essential  for  an 
accurate  solution,  while  allowing  evaluation  of  the 
resulting  integrals  in  a  straightforward  way.  Geometric 
features  which  characterize  a  curve  include  position,  slope, 
curvature,  and  higher  order  derivatives;  but  the  question 
as  to  which  of  these  features  is  essential  for  an  accurate 
solution  has  not  been  adequately  answered  in  the  literature. 

The  approach  which  will  be  used  has  been  considered  by 
Johnson  (Ref  66)  in  a  computer  graphics  context.  This 
approach  is  to  approximate  a  curve  in  a  simple  way  by  using 
a  set  of  standard,  or  primitive,  elements  and  accepting  the 
level  of  error  which  results  from  the  choice  of  the  element. 
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As  noted  by  Johnson,  this  is  in  contrast  to  a  spline 
approach  in  which  geometric  properties  of  a  curve  are 
matched  by  using  as  complex  an  element  as  is  required  to 
accomplish  the  matching. 


The 

simplest 

way 

to  approximate 

a  plane  curve 

of 

moderate 

curvature 

is 

to  use  a  series 

of  straight 

line 

segments 

(Figure 

7). 

Increasing  the 

accuracy  of 

the 

approximation  can 

be 

achieved  by  increasing  the  number 

of 

linear  elements.  The  use  of  higher  order  curves  may  reduce 
the  required  number,  although  at  the  cost  of  introducing 
additional  complexity.  In  an  attempt  to  balance  accuracy 
with  complexity,  several  conic  arc  curves  were  considered. 
In  the  next  section  flat,  parabolic,  circular,  and  elliptic 
arcs  will  be  evaluated  as  to  their  use  as  a  standard 
approximating  element. 

Conic  Arc  Approximation.  Consider  an  arbitrary 
curve  n  =  n(C)  described  in  a  tangent-normal  coordinate 
system  with  the  origin  at  some  point  on  (Figure  8)  so 
that 

n(0)  =  n ' (0)  =  0 

The  problem  is  to  approximate  this  curve  in  the 
region  a  ^  ^  b  using  the  following: 

1.  a  straight  line  segment  given  by  nf(0 

2.  a  parabolic  arc  given  by  np(C) 

3.  a  circular  arc  given  by  n.(£) 
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4.  an  elliptic  arc  given  by  n  (?) 
where  these  curves  are  given  in  the  tangent-normal 
coordinate  system  such  that 

V0)  =  n'h(0)  =0  h  =  f,  p,  c,  e 

An  error  function  E^t?)  can  then  be  defined  to  describe  the 
approximation  in  terms  of  the  error  in  position,  slope,  or 
higher  order  derivatives.  That  is 

e£(?)  =  nn(?)-n£(?) 

where  n  is  the  nth  derivative  of  the  function.  If  these 
curves  are  expanded  in  a  power  series  about  n  =  0  »  and  if 

the  approximating  curves  are  equated  term  by  term  with  the 
actual  curve,  the  errors  of  the  approximation  are  given  by 

Ef(?)  =  0(£2) 

Ep(?)  =  0(C3) 

Ec(5)  =  0(?3) 

E  (?)  =  0(?5) 
e 

Circular  Arc  Approximation.  The  circular  arc  has  been 
chosen  as  the  standard  element  for  several  reasons.  From 
the  analysis  above  it  can  be  seen  that  the  accuracy  of  the 
approximation  increased  as  the  nature  of  the  element  changes 
from  flat  to  elliptic,  and  that  the  errors  resulting  from 


Numerical  Implementation.  The  problem  is  t'  >del  a 
planar  curve  using  a  piecewise  continuous  set  of  circular 
arc  elements.  It  is  assumed  that  the  curve  can  not  be 
described  analytically,  and  will  be  represented  by  a  set  of 
coordinate  points,  input  in  the  case  of  an  airfoil  from  the 
trailing  edge  in  the  order  indicated  in  Fig  9.  The  general 
equation  of  a  circle  is 

x2+y2+Dx+Ey+F  =  0 

where  D,  E,  and  F  are  constants;  where  the  center  of  the 

circle  is  located  at  x0  =  ,  and  yo  -  ^  >  and  where 

the  radius  of  the  circle  is  a=/t>2/4  +  E2/4-F'  •  Given 

three  points  on  the  surface  (Fig  10)  these  three  constants 

can  be  obtained.  Let  P  (x  ,y  )  for  1x1=1, 2, 3  be  three 

m  m  m 

points  on  the  surface  such  that 

x  D+y  E+F  =  -(x  2+y  2 )  m  =  1,2,3 
mm  m  m  ’  * 

This  system  can  be  solved  to  give 

D  =  Ix32+y32-(xi2+y12)3(y2-y3 )-(y !-y3 ) [x3 2+y3 z-(x2 2+y? 2 )1 

W 

E  =  1XX-X3 )[x32+y32-(x22+y22 )]-(x2-x3 )[x32+y32-(xi2+y, 2)1 

W 

F  =  -(x22+y3 2 )-x3D-y3E 
where 

w  =  (xi-x3 )(y2-y3 )-(x2-x3 )(yi-y3 ) 
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Figure  9  -  Surface  Coordinate  Order  of  Input 


Input  Points 


x 


Figure  10  -  Circular  Arc  Panel  Geometric  Definition 


and  from  these  the  circle  radius  and  the  location  of  the 


center  can  be  obtained. 

The  computer  program  which  implements  this  procedure 
sets  an  arbitrary  lower  bound  for  |w|  which  effectively 
determines  how  close  the  three  points  can  be  to  lying  on  a 
straight  line.  While  a  straight  line  can  be  interpreted  as 
a  circle  of  infinite  radius,  the  computer  code  will  not 
accept  this. 

The  effects  of  element  curvature  can  be  studied  by 

passing  the  approximating  circular  arc  through  three 

points  Pi  ,  P2  ,  and  P  (Figure  11)  rather  than 

through  Pi  ,  P2  ,  and  P3  as  described  above.  The  point  P 

is  defined  by  a  parameter  6  given  by 

in 

0  <  8  *  jE.  <  1 

where  l  and  S,  are  defined  in  Figure  11.  As  B  varies 
from  0  to  1,  P  moves  from  P<*  to  P3  ,  and  the  curvature  of 
the  resulting  circular  arc  changes  from  0  to  the  curvature 
of  the  circle  through  the  original  input  points  Pi  ,  p2 
and  P3  .  Note  that  the  approximating  arc  always  goes 
through  the  points  and  P2  .  The  result  is  that  the 

actual  airfoil  is  modeled  as  a  series  of  connected  circular 
arcs.  Having  developed  an  approximation  to  the  surface 

geometry,  the  next  step  in  the  formulation  of  a  panel  method 
is  the  representation  of  the  surface  singularity 
distribution.  This  will  be  discussed  in  the  next  section. 
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Singularity  Representation 


It  has  been  shown  that  source,  doublet,  or  vortex 
singularities  can  be  used  to  model  potential  flow  problems, 
and  that  the  doublet  and  vortex  singularities  are 
equivalent.  The  source  singularity  is  incapable  of 
generating  lift  on  a  body  and  its  exclusive  use  would  be 
unsuitable  for  lifting  cases.  Beyond  this,  however,  there 
is  little  information  available  to  indicate  which 
singularity  is  the  better  one  to  use  for  particular 
applications.  In  terms  of  modeling  the  physical  flow  it  is 
felt  that  the  vortex  singularity  is  more  directly  related  to 
the  actual  flow  since  one  is  trying  to  model  the  viscous 
effects  of  the  boundary  layer  by  using  the  potential  vortex 
sheet  on  the  surface.  The  source  or  doublet  singularity  is 
more  difficult  to  interpret  physically,  and  for  this  reason 
it  is  believed  that  the  vortex  singularity  provides  more 
insight  into  what  is  actually  happening  in  the  flow  near  the 
surface.  Since  part  of  the  purpose  of  this  effort  is  to 
study  the  effect  of  choice  of  singularity  on  the  solution 
the  use  of  both  the  source  and  the  vortex  singularities  will 
be  investigated. 

Series  Expansion.  The  singularity  strength 
distribution  will  be  represented  as  a  series  expansion  of 
the  form 
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(47) 


Q-l 

cr  (  0  )  =  qk  sink(0-01) 

k=0 

while  for  the  airfoil  cases  the  series  used  was  an  expansion 
about  the  panel  midpoint  8„  ,  so  that 

irl 

c(9)  =  ^2  qk  sink(e-9ji)  (48) 

k=0 

Here  Q  is  the  number  of  terms  in  the  series.  In  the 
applications  to  follow  Q  will  be  varied  from  1  to  4  in  order 
to  study  its  effect  on  the  solution.  It  should  be  noted 
that  Q  is  also  the  total  number  of  unknowns  on  a  panel  so 
that  the  total  number  of  unknowns  for  a  problem  which  is 
modeled  by  N  panels  will  be  Q*N .  The  two  forms  for 
singularity  strength  given  above  are  essentially  equivalent, 
although  they  exhibit  certain  differences  in  numerical 
characteristics  which  will  be  discussed  in  later  sections 
where  they  are  applied. 

Continuity  Conditions.  An  important  part  of  this  study 
will  be  to  consider  the  effects  of  singularity  strength 
continuity  on  the  accuracy  of  the  solution.  This  can  be 
done  by  numerically  requiring  continuity  of  the  singularity 
strength  and  its  derivatives  across  panel  junctures.  The 
purpose  of  specifying  continuity  of  derivatives  is  to  obtain 
continuity  of  slope  (or  higher  derivatives)  of  the  function 
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TABLE  IV 

Definition  of  Continuity  Parameter  C 


c 

Singularity  characteristic  at  panel  junctures 

0 

discontinuous 

1 

continuous 

2 

continuous  derivative 

3 

continuous  2nd  derivative 

with  respect  to  arc  length  along  the  airfoil.  It  will  be 
convenient  to  characterize  the  singularity  strength  in  terms 
of  the  degree  of  continuity  which  is  imposed.  For  this 
purpose  a  continuity  class  parameter,  C,  can  be  defined  as 
in  to  Table  IV. 

On  a  circular  arc  of  radius  a  the  arc  length  S  is  given 
by  s  =  a0  so  that 

5no  =  1_  3^o 

Ss11  a11  aen 

Now  if  (  )  refers  to  quantities  on  the  jth  panel,  then 

J 

the  condition  of  continuity  of  the  function  from  panel  j  to 
panel  j  +  1  (that  is  a  class  C  =  1  function)  can  be 

expressed  using  Eq  48  as 


65 


(49) 


aj(eMj+6j)  "  aj+i(9:.iJ+1'  <SJ+1) 

Likewise  the  condition  for  continuity  of  derivatives  can  be 
expressed  as 


i  3n<w v  .  i  alVi(lVi-«.i+i: 


(50) 


Equating  these  expressions  gives  equations  which  enforce 

continuity  of  the  function  across  panel  junctures  through 

the  derivative.  The  singularity  strengths  can  be 

written  using  Eq  48  as 

Q-l 

aj(0M.+<5j)  =  ^2  qk.  sink(6j>  Q  38  1>2>3>  or  4 
J  k=0  j 


+  +  +  -  L  qkJ+1Sll‘k(5j+l) 


Substituting  these  into  Eq  49  and  50  gives  a  matrix  equation 
which  prescribes  continuity  to  class  C  at  panel  junctures  of 
the  form 


Q-l 

Z[cki<v  - 

k=0 


(51) 


Here  the  [C^]  denote  NXN  continuity  coefficient  matrices 
and  the  {Q^}  denote  NX1  column  vectors  whose  elements  are 

p 

the  unknown  q^'s  •  The  derivation  of  the  [C^]  matrices 
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procedure  for  reducing  the  velocity  boundary  condition  to  a 
matrix  equation  will  be  given,  followed  by  a  discussion  of 


the  Kutta  condition  which  will  be  used. 


Velocity  Boundary  Condition, 


To  illustrate  the 


procedure  of  reducing  the  boundary  condition  to  a  matrix 
equation,  consider  the  case  of  a  normal  velocity  boundary 
condition  applied  to  a  surface  on  which  is  placed  a  source 
distribution.  From  Chapter  II  (Eq  27)  the  induced  normal 
velocity  at  a  point  p  on  the  surface  is  given  by 


vn(p> 


_  lim  1  £  .  .  3Ks 

=  P~p  27  J  a(q)  37"  dl 


which  becomes 


vn(p)  -  -  y  /  0(q)  ay 


Now  if  the  surface  is  modeled  with  N  panels,  and  if  the 
source  distribution  on  panel  j  is 

Q-l 

°j(q)  =  2  fk1  (q) 
k=0  J  3 


where  the  f,  are  known  functions  and  the  a  are 

kj  kj 

unknown  constants,  then  the  normal  velocity  at 

the  itla  control  point  is 


N  Q-l 


Vpi> 


2E  P-Pt  2-  J  qkfk  3n  (Ks  ,( 

J-l  J  J  Pi  1J 


panel ^ 


The  integral  over  the  surface  becomes  a  sum  of  integrals 
over  each  panel  which  can  be  computed  analytically.  In 
matrix  form  this  can  be  written 


<v  -  E  CEk]{v 


where  [Rk]  are  called  aerodynamic  influence  coefficients 
and  are  given  by 


?  -  lim  1 

kij  "  P^i  27r 


panel 


3(K_  ) 

fk  — dli 
J  ^  9npi  J 


v»  (Pi) 
n 


is  the  normal  velocity  at  the 


control  point  due  to  the  free  stream,  the  statement  of  zero 
total  normal  velocity  becomes 

Q-l 

E  CRk]  {v  ■  <-v. » 

k=0  n 

The  detailed  equations  for  the  [R^]  are  given  in  Appendix 
B.  The  formulation  if  a  vortex  singularity  is  used  is  the 
same,  and  influence  coefficients  for  this  case  are  also 
given  in  Appendix  B.  As  noted  earlier  the  vortex  case 


requires  an  additional  condition  to  ensure  uniqueness  of  the 
solution . 

Kutta  Condition.  The  requirement  for  a  Kutta  condition 
in  a  lifting  potential  flow  model  is  well  known.  In  Chapter 
II  it  was  noted  that  such  a  condition  was  needed  to  obtain 
a  unique  solution  to  the  problem  of  flow  over  a  body  using  a 
surface  vorticity  distribution.  The  application  of  the 
Kutta  condition  is  an  important  step  in  a  potential  flow 
model  because  it  is  essentially  the  link  between  the  real 
viscous  flow  and  the  potential  model  that  allows  an  accurate 
determination  of  the  lift  on  a  body. 

Theoretically  the  Kutta  condition  requires  a  finite 
flow  velocity  at  a  sharp  trailing  edge.  While  there  are 
many  ways  to  achieve  this  requirement,  three  methods, 
depicted  in  Figure  13,  have  been  used  in  this  study;  a 
specification  of  net  circulation,  a  trailing  edge  bisector 
condition,  and  a  specification  of  zero  vorticity  at  the 
trailing  edge.  The  specification  of  net  circulation  will  be 
used  to  study  the  circular  cylinder  problem,  but  it  is  not 
useful  in  the  study  of  a  general  airfoil  since  the  net 
circulation  is  not  known. 

The  trailing  edge  bisector  condition  (Figure  13b) 
involves  a  specification  of  zero  velocity  at  a 
point  Ax  from  the  airfoil  trailing  edge,  and  in  a 
direction  normal  to  a  line  which  bisects  the  airfoil 


a.  Specification  of  Met  T 


b.  Trailing  Edge  Bisector 


c.  Specification  of  y  at  Trailing  Edge 
Figure  13  -  Kutta  Condition  Specification 


trailing  edge  angle.  The  distance  Ax  ,  which  must  be 
specified  a  priori,  should  be  small,  but  beyond  that  it  is 
arbitrary.  While  this  procedure  has  given  good  results  in 
solving  the  airfoil  problem,  the  arbitrariness  in 
choosing  Ax  is  undesireable . 

The  last  method,  specification  of  zero  vorticity  at  the 
trailing  edge,  has  been  used  to  obtain  the  majority  of  the 
airfoil  results  which  will  be  presented  later  because  it  is 
felt  that  it  is  conceptually  the  most  logical  approach.  It 
is  based  on  the  fact  that,  for  a  surface  distribution  of 
vorticity  with  an  appropriate  boundary  condition  such  that 
the  internal  flow  is  stagnated,  the  surface  vorticity  equals 
r  the  tangential  velocity  on  the  surface.  Since  the 
trailing  edge  should  be  a  stagnation  point  where  the 
tangential  velocity  is  zero,  a  specification  of  zero 
vorticity  at  the  trailing  edge,  both  upper  and  lower 
surface,  is  an  equivalent  Kutta  condition.  A  problem  with 
this  specification  is  that  two  equations  are  required,  one 
for  the  upper  surface  and  one  for  the  lower  surface,  and 
thus  the  complete  problem  is  overspecified  by  one  equation. 

One  way  to  circumvent  this  is  to  use  an  error  parameter 
approach  (Refs  26,  31)  in  which  a  uniform  but  unknown  error 
t  in  normal  velocity  is  assumed  at  each  control  point.  The 
equation  for  zero  normal  velocity  becomes 


t  +  V 
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where  V  is  the  ith  row  of 
ni 

Q-l 

<V  -  E  [V<V 

k=0 

More  generally  this  can  be  written  as 

[M] {T}+{V  }  =  {-Vm  } 
n  <»_ 


and  E  is  a  diagonal  matrix  of  weighting  factors  which 
selects  the  control  points  at  which  the  error  will  be 
applied.  Usually  the  nonzero  elements  of  M  would  be  1. 
Henshaw  and  Bristow  have  had  success  with  this  approach,  but 
it  is  felt  that  the  major  drawback  to  it  is  the  high  degree 
of  arbitrariness  it  introduces  into  the  formulation. 

In  a  second  method,  used  by  Woodward  (Ref  38),  a  source 
of  unknown  strength  is  placed  inside  the  airfoil  to  provide 
the  required  additional  unknown.  Recalling  the  discussion 
of  potential  theory  it  was  found  that  the  problem  of  a 
vortex  distribution  on  a  closed  body  with  an  external  normal 
velocity  boundary  condition  produced  an  ill  posed  problem 
unless  the  net  inflow  through  the  surface  was  zero. 
Theoretically  this  is  the  problem  under  consideration,  but 
in  the  numerical  formulation,  the  net  inflow  condition 
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cannot  be  satisfied  because  boundary  conditions  are  applied 
only  at  discrete  points.  The  internal  point  source  can  be 
thought  of  as  a  method  of  removing  any  excess  fluid  that 
flows  into  the  body  as  a  result  of  imperfect  satisfaction  of 
the  boundary  condition.  The  procedure  is  to  place  a  point 
source  inside  the  body,  and  add  a  term  reflecting  the  effect 
of  the  source  to  the  equation  for  normal  velocity  at  each 
control  point.  Since  the  location  of  the  source  is 
arbitrary,  this  is  a  parameter  whose  effect  on  solution 
accuracy  must  be  studied. 

Numerical  Implementation 

The  basic  elements  of  a  new  panel  method  have  been 
developed  in  the  preceeding  sections.  They  include  the 
choice  of  surface  and  singularity  representations,  the 
selection  of  velocity  and  continuity  boundary  conditions, 
and  if  necessary  the  choice  of  a  Kutta  condition.  In  this 
section  these  elements  will  be  combined  into  a  system  of 
linear  algebraic  equations  which  will  be  solved  using 
standard  methods. 

Matrix  Equation  Formulation.  The  procedure  for 
obtaining  the  solution  to  the  airfoil  problem  using  a  source 
distribution  is  to  model  the  airfoil  with  N  panels.  The 
number  of  terms,  Q,  in  the  singularity  representation  is 
then  chosen,  so  that  the  problem  has  a  total 
of  q»n  unknowns.  The  desired  continuity  class,  C,  of  the 
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singularity  strength  is  then  chosen.  This  formulation  meai-s 
that  with  Q  unknowns  per  panel  and  C  continuity  conditions 
per  panel,  a  total  of  L  =  Q  -  C  velocity  boundary  conditions 
must  be  applied  on  each  panel  in  order  to  have  a 
determinate  linear  algebraic  system.  These  control  points 
will  in  general  be  equally  spaced  on  a  panel,  as  shown  in 
Figure  14. 

The  full  system  of  normal  velocity  boundary  conditions 
and  continuity  conditions  can  be  written 


£ 


k=0 


h  =  1,2, ...L 


Q-l 

J]  [C£]{Qk}  ='o  m  =  0,1, ...C 
k=0 


These  can  be  combined  into  one  Q*N  by  Q*N  system 

[ A] {X}  =  {B} 

For  the  vortex  singularity  case  in  which  an  internal  point 
source  is  added,  and  the  vorticity  is  specified  as  zero  at 
the  trailing  edge,  an  additional  unknown  and  one  equation 
must  be  added  so  that  the  full  system 
is  (Q*N+1)  by  (Q*N+1)  . 

Method  of  Solution.  Although  many  algorithms  are 


75 


I 


il 


Control  Point  Location 


ontrol  Points  per 


Panel 


Control  Point  Spacing 


6 


available  for  solving  systems  of  linear  equations,  no 
special  effort  was  made  in  this  study  to  evaluate  different 
methods.  A  standard  routine  from  the  International 
Mathematical  and  Statistical  Library  (IMSL)(Ref  67)  was  used 
in  all  cases.  This  routine  performs  matrix  inversion  using 
Gaussian  elimination  with  equilibration  and  partial 
pivoting.  It  should  be  noted  that  the  system  developed 
above  exhibits  no  special  characteristics  such  as  bandedness 
or  symmetry  which  would  allow  the  use  of  solvers  designed 
for  such  cases. 

In  this  chapter  the  general  panel  method  approach  to 
solving  potential  flow  problems  has  been  outlined  and  the 
details  of  a  new  panel  method  have  been  presented.  The  new 
method  is  based  on  the  use  of  continuous  circular  arc  panels 
to  model  a  two  dimensional  surface.  A  surface  singularity 
represented  as  a  higher  order  sine  series  expansion  is  then 
distributed  on  the  panels.  This  distribution  is  given  a 
specified  degree  of  continuity,  appropriate  velocity 
boundary  conditions  are  applied,  and  the  problem  is  reduced 
to  a  system  of  linear  algebraic  equations  in  which  the 
unknowns  are  constants  in  the  assumed  singularity 
distribution. 

There  are  several  parameters  in  this  formulation  which 
affect  its  results,  including  type  of  singularity,  number  of 
terms  in  the  series,  the  continuity  class,  number  of  panels, 
curvature  of  the  panels,  and  control  point  location  on  the 


panels.  In  the  next  two  chapters  the  method  will  be  applied 
to  the  problems  of  flow  over  a  circular  cylinder,  and  flow 
over  several  different  airfoils.  The  effects  of  the 
parameters  noted  above  will  be  evaluated,  and  the  results 
will  be  compared  with  those  of  other  two  dimensional 


methods. 


Application  to  the  Circular  Cylinder 


The  purpose  of  this  chapter  is  to  apply  the  present 
method  to  the  problem  of  flow  over  a  circular  cylinder, 
which  has  been  chosen  as  a  test  case  for  a  number  of 
reasons.  First,  the  circle  is  a  simple  shape  for  which  the 
exact  solution  in  terms  of  both  singularity  distribution  and 
surface  velocity  is  easily  computed.  Additionally,  the 
surface  is  free  of  slope  discontinuities  which  will  remove 
the  ambiguities  noted  earlier  which  are  associated  with  a 
surface  singularity  distribution  at  a  corner  and  with  the 
application  of  the  Kutta  Condition.  It  is  realized,  of 
course,  that  the  circle  would  seem  to  be  ideally  suited  for 
a  method  which  uses  circular  arcs  for  panels. 

An  extensive  study  will  also  be  conducted  to  determine 
the  effect  of  a  number  of  parameters  on  the  accuracy  of  the 
solution.  These  include  the  singularity  type,  the  number  of 
terms  in  the  singularity  distribution,  the  number  of  panels, 
continuity,  and  control  point  location.  The  following 
sections  will  discuss  the  exact  solution  with  which  che 
computed  solution  will  be  compared,  the  panel  method 
formulation  of  the  problem,  and  the  results  of  the  parameter 
sensitivity  studies. 
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Exact  Solution 


Consider  the  classic  problem  of  uniform  flow  over  a 
circular  cylinder  as  shown  in  Figure  15.  Although  this 
problem  can  be  solved  in  several  ways,  the  exact  solution 
will  be  developed  in  terms  of  source  and  vorticity  surface 
singularity  distributions.  This  will  provide  an 
introduction  to  the  use  of  this  method  to  obtain  approximate 
solutions  to  more  complicated  problems. 

Source  Distribution.  Assume  there  exists  a  surface 
source  distribution,  a(9)  ,  on  the  cylinder  shown  in 
Figure  15,  and  apply  a  zero  normal  velocity  condition  to 
this  problem.  Let  P  go  to  p  on  the  surface  r=l  to  obtain 
from  Eq  28,  the  normal  velocity  component  induced  by  the 
source  sheet  as 


vr(l,B) 


mm  " 

=  +  1?  /  a(9  o  )d0  o 


(52) 


The  boundary  condition  is 


vr(l,0)+Voo  -  vr(l,0)+cose  =  0 


so  that  Eq  52  becomes 


*•  II 

<?(9)  +  ^  J  a(9o)d0o  “  -2  cos  9 


(53) 


This  is  a  Fredholm  Equation  of  the  Second  Kind  with 
parameter  X  ■  .  Since  X  is  not  an  eigenvalue  Eq  53  has 
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Circular 

Cylinder 


Figure  15  -  Flow  Over  a  Circular  Cylinder 


a  unique  nontrivial  solution, 
integrating  Eq  53  to  get 


It  can  be  obtained  by 


Letting 

2tt 

%  -  h  f  0<e°>d6« 

0 

where  q  is  the  total  source  strength 
s 

Eq  54  becomes 

2ttQs+Qs  =  0 

Therefore 

Qs  =  0 

and 

a(9)  =  -2  cos  9 

The  induced  normal  and  tangential  velocities  on  the  surface 
are 

vr(l,9)  =  -cos9 

and 

v  (1,9)  =  -sin9 

and  the  total  tangential  surface  velocity  is 
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V.  =  v.  +  V  =  -2  sin  0 
*  4  t 


(56) 


The  natural  result  of  this  formulation  is  that  the  total 
source  strength  is  zero,  as  should  be  expected  both  in  this 
exact  solution  and  in  the  subsequent  approximate  solutions. 
The  deviation  from  zero  of  the  total  source  strength  can  be 
used  as  a  measure  of  the  accuracy  of  the  approximate 
solution. 

As  was  discussed  in  Chapter  II,  the  flow  in  the 
interior  region  of  a  closed  body  is  independent  of  that  in 


the  exterior  region.  For  the  choice  of 


source 


distribution  on  the  cylinder  with  a  zero  normal  velocity 
boundary  condition,  the  flow  pattern  in  the  interior  of  the 
cylinder  is  shown  in  Figure  16a. 

Vorticity  Distribution.  This  problem  can  also  be 
solved  using  a  surface  vorticity  distribution  and  a  zero 
internal  tangential  velocity  boundary  condition.  In  this 
case  the  perturbation  tangential  velocity  on  the  interior 


surface  of  the  cylinder  due  to 
distribution,  Y(9)  is 


vorticity 


vt(l,0) 


-Y(0' 


A  I 

hi 


Y( 6  o )d0  o 


(57) 


Applying  the  boundary  condition 


vt  +  Vt  =  v^.  -  sin  0  =  0 
00 
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Eq  57  becomes 


td  ft 

y(9)-  J  y( 0  o )d6  o  =  -2  sin  6 


(58) 


Recall  that  this  is  a  Fredholm  Integral  Equation  of  the 
Second  Kind  with  parameter  X  =  ^  .  Since  X  is  an 
eigenvalue,  a  solution  to  Eq  58  exists  only  if 


mm  II 

/  u( 9  0  )f ( 


0  o )d0  o  —  0 


(59) 


when  f ( 0  o )  =  -2  sin  0O 

and  where  u(0o)  is  a  solution  to 


cm  n 

u(0)  ^  j  u(0  o  )d0o  =  0 


(60) 


Clearly  u(0o)  *  c. 


(where  c .  is  any  constant  is  a 
J 


solution  to  Eq  60,  and  thus  is  an  eigenfunction  of  Eq  58. 

Given  any  c.  Eq  59  holds;  therefore,  the  general  solution 
J 

to  Eq  58  is 


Y ( 9 )  =  Y„( 0 )+  2  c, 

where  y  (0)  is  a  particular  solution  and  E  c.  is  a 

P  j  3 

sum  over  all  linearly  independent  eigenfunctions.  But  since 

each  Cj  is  a  constant  only  one  will  be  linearly 

independent;  therefore  z  c=D  where  D  is  any  constant. 

j  3 

By  inspection,  y_(9)  =  -2  sin  0  ,  so  that 


inside  the  cylinder  (see  Figure  16b)  compared  with  the  non 
zero  flow  produced  by  the  source  distribution. 


Combined  Source /Vorticity  Distribution.  The  preceeding 
two  approaches  can  be  combined  by  assuming  source  and 
vorticity  distributions  of  the  form 

a(9)  =  -A  cos  9  (64c) 

y( 9 )  =  -B  sin  9  +  i—  (64b) 

Z7T 

where  A  and  B  are  constants  to  be  determined.  The  total 
external  normal  velocity  is  then 


Vn(l,0)  »  (1-  ^2)  cos  9  (65) 

and  a  zero  normal  velocity  boundary  condition  requires  that 

A  +  B  -  2  (66) 

The  total  tangential  velocity  is  then,  using  Eq  66, 

Vt(l,9)  =  -sin  9  -  sin  9  +  ^  =  -2sin9+  ^  (67) 

The  velocity  given  by  Eq  67  is  the  same  for  any  values  of  A 
and  B  as  long  as  Eq  66  is  satisfied. 

For  example,  choosing  A  =  B  =  1  is  equivalent  to 
setting  the  source  strength  equal  to  the  normal  component 


Panel  Method  Solution 


The  method  developed  in  Chapter  III  will  now  be  applied 
to  the  circular  cylinder  problem.  The  circle  will  be 
divided  into  N  panels  of  equal  arc  length  with  panel  number 
1  centered  on  the  trailing  edge  stagnation  point,  as  shown 
in  Figure  17.  While  Eq  48,  termed  the  element  centered 
formulation,  is  the  preferred  singularity  distribution  for 
the  case  of  a  general  body  ,  when  it  is  used  for  the  case  of 
a  circle  with  equally  spaced  panels  with  control  points  at 
panel  centers  (the  points  about  which  the  distribution  is 
expanded)  some  elements  of  the  velocity  influence 
coefficients  become  zero,  producing  a  singular  matrix.  For 
this  reason  Eq  47,  termed  the  element  non-centered 
formulation,  will  be  used  to  represent  the  singularity 
distribution  for  flow  over  a  circle. 

An  advantage  of  using  this  series  is  that  the 
continuity  matrices  become  diagonal.  A  disadvantage  is  that 
the  results  are  not  completely  symmetric.  The  degree  of 
symmetry  increases  as  the  overall  accuracy  of  the  solution 
is  increased  by  varying  other  parameters.  The  results  do 
exhibit  a  polar  symmetry  about  the  origin.  That  is,  the 
results  on  a  ray  connecting  two  points  on  the  circle  and 
passing  through  the  origin  are  identical.  Note  that  a 
solution  which  assumes  symmetry  has  not  been  developed  so 
that  the  method  may  be  applied  to  asymmetric  airfoils. 
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Results 

The  results  of  the  application  of  the  method  to  the 
circular  cylinder  problem  will  be  presented  in  terms  of 
velocity  error  plots  since  the  exact  velocity  on  the 
cylinder  can  be  computed.  First,  the  effects  of  the  various 
parameters  will  be  compared  using  a  series  of  maximum 
absolute  velocity  error  plots.  This  will  be  followed  by  a 
consideration  of  the  local  error  distribution  on  the 
surface,  and  finally  a  discussion  of  the  sensitivity  of  the 
solutions  to  control  point  location  will  be  presented. 
Although  no  special  attention  was  given  to  the  question  of 
computational  efficiency,  a  limited  assessment  of  the  effect 
of  the  higher  order  method  on  efficiency  will  be  made.  In 
general,  each  solution  presented  required  no  more  than 
several  seconds  of  computer  time  on  a  CDC  6600/CYBER  74 
computer. 

Global  Error.  Figures  18  to  21  show  the  effects  on 
accuracy  due  to  panel  size,  or  number  of  elements  (N), 
continuity  (C),  and  number  of  terms  (Q)  in  the  singularity 
distributions  for  various  choices  of  singularity  and 
boundary  condition.  Figure  18  also  shows  lines  of  constant 
computer  time  which  will  be  discussed  later.  These  charts 
show  the  maximum  absolute  value  error  in  surface  normal  or 
tangential  velocity  versus  N  for  various  combinations  of  Q 
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MAXIMUM  ERROR  IN  VN 
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Figure  18.  Maximum  Absolute  Normal  Velocity  Error  on  a 
Circle,  Source  Distribution,  Normal  Velocity 
Boundary  Condition. 
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Figure  19. 


Maximum  Absolute  Tangential  Velocity  Error  on 
a  Circle,  Source  Distribution,  Normal  Velocity 
Boundary  Condition. 
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Figure  20.  Maximum  Absolute  Normal  Velocity  Error  on  a 
Circle,  Vortex  Distribution,  Normal  Velocity 
Boundary  Condition. 
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Figure  21.  Maximum  Absolute  Tangential  Velocity  Error  on 
a  Circle,  Vortex  Distribution,  Normal  Velocity 
Boundary  Condition. 


and  C.  The  ranges  on  these  parameters  are 


4  <  N  £  48 

1  -  Q  -  4 

0  <  C  <  3 

L  +  Q  -  C  -  Control  Points 

Per  Panel 

From  the  solution  obtained  for  each  case,  velocities  on  the 
circle  were  computed  at  120  equally  spaced  points  around  the 
circle,  and  from  these  the  largest  absolute  value  errors 
were  determined  according  to  Eq  69. 


Eqc  =  MAX 


V  (0 . )  -  V  (8 . ) 
comp^  i'  exv 


(69) 


1,2 


.120 


This  approach  is  a  simple  way  of  comparing  the  effect  of  the 
above  parameters  on  the  relative  accuracy  of  the  computed 
solutions  for  various  choices  of  the  parameters,  and  it 
allows  the  effects  of  these  parameters  on  the  solution  to  be 
studied.  The  values  for  maximum  error  are  not  to  be 
interpreted,  however,  as  the  largest  error  anywhere  on  the 
surface  for  a  particular  solution.  Since  the  velocities  are 
singular  at  panel  endpoints  the  error  there  can  be  made 
arbitrarily  large  by  computing  velocities  at  points  closer 
and  closer  to  the  panel  endpoints  (at  least  for  the 
discontinuous  cases).  Another  fact  to  note  concerning  the 
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above  figures  is  that  points  are  missing  for  certain 
parameter  combinations  because  the  influence  coefficient 
matrices  were  singular.  Examination  of  the  matrices  in 
question  revealed  that  this  phenomena  is  a  numerical  result 
of  the  symmetry  of  the  circle  problem. 

Figures  18  and  19  show  the  maximum  absolute  errors  in 
normal  and  tangential  velocity  respectively  for  a  source 
distribution  on  the  circle  using  a  normal  velocity  boundary 
condition.  Figures  20  and  21  show  the  same  maximum  absolute 
errors  for  a  vorticity  distribution  for  which  the  total 
circulation  was  specified  as  zero.  Results  for  the  case  of 
non  zero  circulation  were  nearly  identical. 

The  symmetry  of  the  problem  resulted  in  three 
interesting  effects.  First,  the  total  source  strength  is 
identically  zero  for  all  cases  as  it  should  be  for  an  exact 
solution.  Second,  total  circulation  in  the  vorticity  case 
can  be  specified  without  adding  an  additional  equation  to 
the  system.  The  reason  for  this  is  that  the  expression  for 
net  circulation  is  embedded  in  the  left  hand  side  of  the 
problem,  and  can  be  conveniently  extracted  and  fixed. 
Although  this  would  not  in  general  be  an  acceptable  method 
of  specifying  circulation,  the  symmetry  of  the  problem 
ensures  that  this  technique  will  be  successful  for  the 
circular  cylinder  case.  Third,  while  the  case  of  a  vortex 
only  distribution  with  normal  velocity  boundary  conditions 
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difference,  and  the  rate  of  accuracy  improvement  both 
increase  up  to  the  Q=3  cases.  However  as  Q  goes  from  Q=3  to 
Q=s4  the  effects  of  continuity  seem  to  play  a  more  important 
role  than  in  the  source  case.  For  the  Q=2  case  the 
continuous  solution  is  much  better  than  the  discontinuous 
solution,  in  contrast  to  the  source  case  where  the  two  were 
very  close,  with  the  discontinuous  case  being  slightly 
better.  Considering  the  Q=3  cases,  it  is  found  that  both 
the  C=0  and  C=2  cases  are  actually  better  than  the  Q=4,  C=0 
case.  In  fact,  if  E..,-,  denotes  the  error  for  the  case  Q 
and  C,  the  ylative  level  of  error  for  the  normal  velocity 
is  seen  to  be 


J41 


<  E/io  <  <  Eoo  <  E 


43 


30 


32 


■"40 


while  that  for  the  tangential  velocity  is  seen  to  be 


J41 


<  E 


43 


<  E 


32 


<  E30  <  E 


40 


For  a  given  value  of  Q,  continuity  is  important  for  the 
tangential  velocity  error,  but  additional  degrees  of 
continuity  beyond  C=1  are  not  required.  Also,  the  effect  of 
increasing  Q  from  3  to  4  is  not  clear  in  that,  contrary  to 
what  might  be  expected,  the  benefit  of  the  additional  term 
seems  to  depend  on  the  particular  continuity  requirement 
that  is  imposed. 

It  has  been  noted  that  a  discontinuous  distribution  of 
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•  Panel  End  Point 
X  Control  Point 


normal  direction 
at  control  point 


sources  gives  better  results  than  does  a  continuous 
distribution,  while  the  opposite  is  true  for  the  vortex 
case.  A  possible  explanation  for  this  lies  in  the 
characteristics  of  the  singularities  themselves  coupled  with 
the  normal  velocity  type  boundary  condition  which  was  used 
in  the  above  cases. 

For  the  case  of  a  vortex  distribution,  a  discontinuity 
in  strength  at  a  panel  juncture  will  act  like  a  line  vortex, 
the  effect  of  which  will  be  felt  mainly  by  control  points 
near  the  panel  juncture,  as  opposed  to  control  points  which 
are  far  away  from  it.  The  velocity  induced  by  this  line 
vortex  at  these  nearby  control  points  will  usually  have  a 
significant  component  normal  to  the  panel  because  most 
panels  are  not  highly  curved.  This  can  be  seen 
qualitatively  in  the  sketch  above.  This  normal  component  of 


normal  direction 
at  control  point 


velocity  induced  by 
source  discontinuity 


velocity  which  would  not  be  there  if  the  singularity 
distribution  were  continuous,  is  very  effectively  cancelled 
by  the  normal  velocity  boundary  condition  which  was  used, 
but  the  vortex  strength  solution  thus  obtained  is  not  what 
it  would  be  if  the  vortex  distribution  was  continuous.  Thus 
it  seems  reasonable  that  the  vortex  case  would  be  sensitive 
to  whether  or  not  a  continuous  distribution  was  used. 

It  is  also  reasonable  to  expect  that  the  source 
distribution  might  be  less  sensitive  to  imposition  of 
continuity  for  similar  reasons.  In  the  source  case  the 
velocity  component  induced  by  a  discontinuity  at  a  panel 
juncture  at  a  nearby  control  point  would  be  small  if  the 
panel  was  not  highly  curved  (See  the  sketch  above).  Thus 


this  normal  component  of  velocity  would  not  have  a  large 


effect  on  the  source  strength  solution  obtained  by  applying 
a  normal  velocity  boundary  condition. 

Hess  (Ref  12)  has  developed  a  criterion  for 
mathematical  consistency  of  a  panel  method.  This  criterion 
is  that  the  singularity  distribution  should  be  of  an  order 
one  degree  lower  than  the  order  of  the  surface  element.  He 
notes,  however  that  others  (Ref  25)  have  violated  this  rule 
and  have  obtained  good  results.  The  present  results  for  the 
cylinder  given  above  are  also  in  violation  of  this  rule, 
because,  although  the  circular  arc  is  a  quadratic  element, 
accuracy  improvements  were  obtained  for  the  Q=2,3  and  4 
cases.  It  is  felt  that  the  errors  introduced  during  the 
numerical  implementation  of  the  present  method,  or  other 
methods,  probably  overshadow  the  mathematical  argument  for 
consistency. 

Figure  22  shows  the  maximum  error  in  normal  velocity 
for  a  source  versus  a  vortex  solution  for  several 
combinations  of  Q  and  C.  The  element  number  ranges  from  4 
to  12.  The  vortex  results  were  obtained  using  a  tangential 
velocity  boundary  condition.  The  maximum  error  in  the 
vortex  solution  is  seen  to  be  consistently  less  than  that 
for  the  source  solution  for  the  same  number  of  elements. 

Efficiency 

Superimposed  on  Figure  18  are  lines  of  constant 


computer  time  for  a  CDC  6600/CYBER  74  computer.  These  times 


hohh3  wwixvw 


0+Dt> 


are  useful  only  for  comparing  the  effects  of  the  parameters 
shown  on  computational  efficiency  for  the  present  method. 
It  can  be  seen  that  to  obtain  a  given  level  of  accuracy,  for 
example  an  error  of  10E-3,  one  needs  roughly  1.5  seconds  of 
computer  time  for  a  3  term  series,  roughly  2.5  seconds  for  a 
2  term  series,  and  some  much  larger  amount  of  time  for  a  1 
term  series.  This  means  that,  although  for  a  given  number 
of  elements  computer  time  increases  with  Q,  the  level  of 
error  is  decreasing  at  a  faster  rate  than  Q  is  increasing. 
Thus,  for  the  present  method  with  the  ranges  of  Q  and  N 
shown,  the  higher  order  singularity  distribution  is  more 
efficient  than  a  lower  order  distribution. 

Local  Error.  The  previous  discussion  dealt  with  a 
measure  of  what  might  be  called  global  error.  Now  consider 
a  local  error  by  looking  at  the  actual  error  distribution  on 
the  circle  surface.  Figures  23  through  28  show  the  normal 
and  tangential  velocity  errors  on  the  surface,  computed  at 
120  equally  spaced  points,  as  a  function  of  angle  measured 
counterclockwise  from  the  trailing  edge  stagnation  point. 
Results  for  the  zero  circulation  case  are  given  only  for  the 
upper  half  circle  because  of  flow  symmetry.  The  errors 
have  the  form 


00 


The  error  in  normal  velocity  is  a  measure  of  the  leakage 
through  the  surface.  Since  Vnex  *  0  »  Vner  equals 
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Velocity  Error  for 
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Eigure  25b.  Effect  of  C  on  Normal  Velocity  Error  for  a  Circle,  QC-30  and  QC-32 
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Figure  28b.  Effect  of  C  on  Tangential  Velocity  Error  for  a  Circle,  QC=30  and  QC-22 


the  computed  value  Vn  ,  and  it  can  be  shown 
c  comp 

that  2Vn  =  a  -a  where  a  is  the  computed 
comp  ex  comp 

source  strength  and  agx  is  the  exact  source  strength. 
Thus  the  error  in  normal  velocity  is  also  a  direct  measure 
of  the  error  of  the  computed  source  distribution.  All  of 
these  results  are  for  source  distributions. 

Figures  23a  and  23b  show  the  effect  of  Q  on  the  normal 
velocity  error  distribution  for  several  8  element  cases. 
These  are  all  C=0  cases,  and  this  is  reflected  in  the 
discontinuous  nature  of  the  error  distributions.  Note  the 
difference  in  scale  between  the  Q=1  case,  and  the  Q=2  and 
Q=3  cases.  Also  note  the  reduction  in  magnitude  and  the 
general  flattening  of  the  curves  on  a  panel  as  Q  increases. 
This  is  a  result  of  the  fact  that  more  control  points  on  a 
panel  provide  better  control  of  the  normal  flow  through  the 
panel . 

Figure  24  shows  the  effect  of  N  on  the  normal  velocity 
error  distribution  for  the  QC=32  case.  Although  the  overall 
level  of  error  for  either  case  is  small,  the  effect  of 
doubling  the  number  of  elements  is  dramatic.  Although  these 
curves  seem  to  exhibit  an  oscillatory  nature,  this  is  to  be 
expected  since  the  boundary  conditions  are  satisfied  only  at 
discrete  points.  In  between  these  discrete  points  the 
solution  effectively  over-  and  under-shoots  the  correct 
solution.  The  magnitude  of  the  apparent  oscillations 
decreases  as  N  and  Q  increase. 


Figures  25a  and  25b  show  the  effect  of  continuity  on 
the  normal  velocity  error  distribution  on  a  16  element 
circle  for  the  Q=2  and  Q=3  cases.  Although  the  continuity 
of  the  singularity  distributions  is  reflected  in  the 


smoothness 

of 

the  error 

curves  (in 

these 

and 

previous 

figures) 

the 

level  of 

error  does 

not 

seem 

to  be 

significantly  affected  by  the  degree  of  continuity.  This  is 
consistent  with  previous  results  for  the  maximum  absolute 
velocity  errors. 

Figures  26  through  28  present  similar  figures  for  the 
distribution  of  tangential  velocity  error.  The  general 
improvement  in  accuracy  as  N  and  Q  increase  is  apparent ,  as 
is  the  general  flattening  of  the  error  distributions  as  Q 
increases. 

Parameter  Study.  An  extensive  study  using  a  source 
distribution  with  N**16  was  made  to  determine  the  sensitivity 
of  the  solution  to  control  point  location  on  a  panel.  All 
the  results  to  this  point  have  been  for  control  points  which 
were  equally  spaced  on  a  panel.  For  cases  which  require  one 
control  point  per  panel  the  location  of  the  control  point 
was  varied  between  20%  and  80%  of  a  panel's  arc  length.  For 
cases  which  require  more  than  one  control  point  per  panel 
all  but  one  were  equally  spaced  and  fixed  on  a  panel,  while 
one  was  allowed  to  vary  between  the  aforementioned  limits. 
In  one  case  which  required  two  control  points  per  panel  both 
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TABLE  V 

Fixed  Control  Point  Location 
(fraction  of  panel  arc  length) 


Control 
Point  No. 

40 

41 

30 

20 

1 

.25 

.33 

.33 

.50 

2 

Ln 

O 

.67 

.67 

3 

.75 

were  allowed  to  vary. 

Figures  29a  and  29b  show  the  effect  of  control  point 
location  on  maximum  absolute  normal  and  tangential  velocity 
error.  The  format  of  these  figures  is  the  same  as  in 
Figures  18  through  21.  Table  V  shows  the  locations  of  the 
fixed  control  points  for  the  multiple  control  point  cases. 
In  general  these  figures  show  a  relatively  small  effect  of 
control  point  location  on  the  overall  level  of  error  in  an 
absolute  sense,  with  similar  results  for  both  normal  and 
tangential  velocity  errors.  A  general  observation  that  can 
be  made  is  that  the  C=0  curves  are  all  concave  upward,  while 
the  continuous  cases  are  all  concave  downward,  except  for 
the  Q032  case.  This  would  indicate  that  for  C=0  cases 
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Figure  29a.  Effect  of  Control  Point  Location  on  Maximum 
Absolute  Velocity  Error  for  Circle,  Source 
Distribution  with  Different  QC  Combinations, 
and  with  N*16,  Normal  Velocity  Error. 
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Figure  29b.  Effect  of  Control  Point  Location  on  Maximum 
Absolute  Velocity  Error  for  Circle,  Source 
Distribution  with  Different  QC  Combinations, 
and  with  N=*  16 ,  Tangential  Velocity  Error. 
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control  points  nearer  the  panel  center  would  be  more 
effective,  while  the  reverse  would  be  true  for  the  C  not 
equal  to  zero  cases. 


Figure  30  shows  more  detailed  comparisons  for  the  N=16, 
QO20  case.  Since  this  case  requires  two  control  points  per 
panel,  one  was  fixed,  successively,  at  10%,  25%,  and  50%  of 
panel  arc  length,  while  the  other  was  varied  between  20%  and 
90%  of  the  panel  arc  length.  When  the  fixed  control  points 
are  at  10%  and  25%  of  the  panel  arc  length,  the  smallest 
errors  occur  when  the  free  control  point  is  near  80%,  and 
the  minimum  for  the  case  of  the  fixed  control  point  at  50% 
is  when  the  free  control  point  is  also  near  50%.  This 
indicates  that  a  symmetric  placement  of  control  points  on  a 
panel  gives  the  smallest  level  of  maximum  error. 

Figure  31  shows  results  for  the  only  case  is  which  more 
than  one  control  point  was  moved.  In  this  case,  where 
QC-20,  the  control  points  were  initially  placed  at  20%  and 
80%  of  the  panel  arc  length,  and  then  they  were  moved  closer 
together  at  the  same  rate.  The  abscissa,  S,  represents  a 
fraction  of  panel  arc  length  with  one  control  point  at  S  and 
the  second  one  at  1-S.  These  curves  have  minimums  in 
roughly  the  20%  to  30%  range. 

The  results  shown  in  Figures  29a  and  29b  are  absolute 
velocity  errors.  Figures  32a  and  32b  (1  control  point  case) 
and  32c,  and  32d  (more  than  1  control  point  case)  reformat 
this  data  by  normalizing  each  curve  by  its  own  maximum  value 


Figure  30.  Comparison  of  Effect  of  Control  Point  Location 
on  Maximum  Absolute  Normal  and  Tangential  Vel¬ 
ocity  Error  for  a  Circle. 
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Figure  31.  Effect  of 
Velocity 
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Figure  32a.  Effect  of  Control  Point  Location  on  Normalized 
Maximum  Velocity  Error  for  a  Circle,  Source 
Distribution,  N»16,  Normal  Velocity,  1  Control 
Point  Cases. 
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Figure  32b. 


Effect  of  Control  Point  Location  on  Normalized 
Maximum  Velocity  Error  for  a  Circle,  Source 
Distribution,  N=16,  Tangential  Velocity, 

1  Control  Point  Cases. 
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MAXIMUM  ERROR  IN  VN 


125 


SOURCE  DISTRIBUTION 
N  *  16 


□  =QC-20 
a  =QO30 
A=QC=40 
4  =QC=41 


CP  LOCATION  (FRACTION  OF  ARC  LENGTH) 


Figure  32d.  Effect  of  Control  Point  Location  on  Normalized 
Maximum  Velocity  Error  for  a  Circle,  Source 
Distribution,  N=16,  Tangential  Velocity, 
Remaining  Cases. 
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TABLE  VI 

Maximum  Absolute  Normal  and  Tangential  Velocities 


QC 

VN  MAX 

VT  MAX 

10 

280.70955 

225.44110 

20 

31.58174 

21.06721 

21 

17.45405 

10.39468 

30 

.91676 

.44371 

32 

.63308 

.49764 

40 

.20180 

.10712 

41 

.14625 

.08988 

43 

.24076 

.16825 

so  that  the  maximum  value  on  each  curve  is  1.  The  values 
used  to  normalize  each  curve  are  given  in  Table  VI. 
Although  this  presentation  accentuates  the  effects  of 
control  point  location,  previous  comments  remain  valid. 

In  general,  the  variation  in  tangential  velocity  error 
due  to  control  point  location  was  under  40%,  and  was  often 
under  25%.  The  method  was  judged  to  be  not  critically 
sensitive  to  control  point  location  since  shifting  this 
location  did  not  make  order  of  magnitude  changes  in 
accuracy.  The  study  does  not,  however,  indicate 
overwhelming  evidence  which  would  support  one  control  point 
location  over  another  as  a  general  rule.  Thus  the  choice  of 

equally  spaced  control  point  locations  is  a  rational  and 
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acceptable  choice,  and  will  be  used  in  the  airfoil 
applications  which  will  follow. 

Combined  Distribution.  Results  presented  up  to  this 
point  have  been  for  source  distributions  or  for  vortex 
distributions.  Intuitively  one  might  expect  some  advantage 
to  be  gained  by  a  combined  source-vortex  or  source-doublet 
solution  method.  That  this  approach  seems  fundamentally 
sound  can  be  argued  as  follows:  source  singularities  are 
more  effective  near  stagnation  points.  In  the  forward 
stagnation  region  the  free  stream  must  be  countered  by  a 
strong  efflux  and  near  the  rear  stagnation  point  the  flow 
must  be  drawn  in  by  a  strong  influx.  On  the  other  hand, 
vorticity  or  doublet  singularities  are  more  effective  in 
generating  and  controlling  surface  tangential  velocities, 
and  thus  should  dominate  on  those  parts  of  the  body  where 
tangential  velocities  are  large.  This  precise  behavior  is 
demonstrated  by  the  source  only  and  vorticity  only  exact 
solut ions. 

Bristow  (Ref  29)  has  found  that  a  hybrid  method  based 
on  Green’s  third  identity  and  employing  higher  order  curved 
panels  is  both  accurate  and  numerically  stable.  In  this 
formulation,  Green's  identity  specifies  that  the  source 
strength  density,  a  ,  be  equated  to  the  surface  normal 
perturbation  velocity  component  .  Then 


A  similar  hybrid  version  of  the  present  method  was 
developed  by  superimposing  the  source  only  and  vorticity 
only  panel  methods.  This  approach  is  tantamount  to 
splitting  the  free  stream  velocity  into  two  equal  parts,  one 
each  for  solving  source  only  and  vortex  only  problems.  This 
approach  is  identical  to  the  A=B=1  version  of  the 
superimposed  exact  solution  discussed  earlier. 

Solutions  obtained  in  this  way  were  substantially  the 
same  as  the  vorticity  panel  results  and  offered  no  apparent 
advantage  since  numerical  instabilities  were  absent  in  both 
methods.  Figures  33  and  34  show  the  effects  of  N,  Q,  and  C 
on  maximum  absolute  velocity  errors  for  a  source/ vortex 
combined  distribution.  An  attempt  to  deviate  from  the 
Green's  identity  specification  of  a  was  made  by  first 
solving  for  by  the  source  panel  method  using  the  full  free 
stream  velocity  and  then  solving  for  y  under  the  influence 
of  the  source  distribution.  For  both  cyclic  and  acyclic 
problems  this  led  to  =  constant  for  Q=l,  and  numerical 
difficulties  for  Q=2. 

Based  on  these  results  for  the  use  of  circular  arc 
panels  to  model  flow  past  a  circular  cylinder,  with  and 
without  circulation,  the  computational  evidence  indicated  no 
advantage  of  the  hybrid  method,  at  least  when  satisfying 
Neumann  type  boundary  conditions.  Also,  from  the  standpoint 
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Figure  33.  Maximum  Absolute  Error  in  Normal  Velocity  for  a 
Circle,  Source/Vortex  Combination  with  7=0  . 
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Figure  34.  Maximum  Absolute  Error  in  Tangential  Velocity 
for  a  Circle,  Source/Vortex  Combination  with 
r»o  . 
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of  comparing  and  the  physical  vorticity  distribution,  the 
hybrid  method  is  inferior  to  a  vorticity  only  solution. 

Summary.  Several  conclusions  can  be  drawn  from  this 
study  of  the  circular  cylinder  problem: 

1.  The  accuracies  obtainable  from  source  only, 
vorticity  only,  and  combined  source/ vortex  methods  are 
roughly  equivalent.  The  vorticity  method  appears  to  be 
superior  for  three  reasons:  its  applicability  to  flows  with 
lift,  its  more  accurate  results  for  tangential  velocity  and 
thus  surface  pressure,  and  its  more  accurate  modeling  of 
physical  vorticity. 

2.  Dramatic  reductions  in  velocity  errors  are  achieved 
by  increasing  Q,  the  number  of  terms  in  the  series 
representation  of  singularity  strength,  through  Q=3.  A 
linear  distribution  (Q=2)  may,  however,  represent  the  best 
compromise  between  simplicity  and  accuracy. 

3.  Accuracy  improvements  were  achieved  by  increasing  N 
(decreasing  panel  size),  panel  junctures,  +  he  surface 
velocity  distributions  were  smooth  with  continuity  and 
discontinuous  without  it. 

5.  The  method  is  not  critically  sensitive  to  control 
point  location  provided  control  points  are  not  located  too 
near  panel  edges. 
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V.  Application  to  Airfoils 


The  purpose  of  this  chapter  is  to  apply  the  method 
developed  previously  to  several  types  of  airfoils  in  order 
to  assess  its  performance  in  these  cases.  The  airfoils 
which  have  been  studied  are  a  symmetric  Joukowski  airfoil, 
an  NACA  0024  airfoil,  a  thin  symmetrical  airfoil,  and  two 
types  of  Karman-Tref ftz  airfoil.  Initially  a  source 
distribution  was  used  to  compute  the  potential  flow  over  the 
Joukowski  airfoil  to  indicate  whether  or  not  the  kind  of 
results  obtained  in  Chapter  IV  for  the  circle  could  also  be 
obtained  for  an  airfoil.  Since  these  preliminary  efforts 
were  promising,  a  source  distribution  was  then  used  to 
compute  the  flow  over  an  NACA  0024  airfoil.  This  was  used 
because  it  is  more  representative  of  a  real  airfoil  section 
(i.e.  it  has  a  non  zero  trailing  edge  angle  compared  to  the 
Joukowski  airfoil's  cusped  trailing  edge),  and  other 
computational  and  experimental  results  were  available  for 
comparison.  A  vortex  distribution  was  then  used  to  compute 
the  thin  symmetrical  airfoil. 

These  studies  led  to  the  application  of  the  method 
using  source  and  vortex  distributions  to  both  a  symmetric 
and  a  cambered  Karman-Tref ftz  airfoil.  The  bulk  of  the 
effort  was  concentrated  on  these  airfoils  because  they  have 
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a  non  zero  trailing  edge  angle,  and  because  exact  solutions 
are  readily  available  for  comparison.  The  unsuccessful 
results  using  certain  parameter  combinations  in  computing 
these  airfoil  cases  led  to  the  selection  of  a  baseline 
method  which  was  then  used  in  an  extensive  parameter  study 
and  error  analysis.  This  study  was  made  for  both  types  of 
Karman-Tref ftz  airfoils  at  various  angles  of  attack.  It 
included  the  effects  of  panel  size,  number  of  terms  in  the 
singularity  distribution,  panel  distribution  characteris¬ 
tics,  and  point  source  location.  An  analysis  of  the  error 
introduced  into  the  solution  by  the  error  in  surface  slope 
and  position  at  control  points  was  also  accomplished. 

In  the  succeeding  sections  of  this  chapter  results  of 
the  application  of  the  method  to  a  Joukowski  airfoil  will  be 
presented,  followed  by  results  for  the  NACA  0024  airfoil, 
and  for  a  thin  symmetrical  airfoil.  Results  for  Karman- 
Tref  ftz  airfoils,  including  detailed  parameter  studies,  will 
then  be  presented.  Finally,  conclusions  regarding  the  value 
of  the  method  for  airfoil  applications  will  be  discussed. 

Joukowski  Airfoil 

As  an  initial  test  of  the  method  when  applied  to  an 
airfoil,  a  13%  thick,  symmetric  Joukowski  airfoil  was 
modeled  using  a  source  distribution.  Cases  were  computed 
for  several  combinations  of  N  and  Q,  and  the  effect  of  panel 
curvature  was  also  investigated. 
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TABLE  VII 


arrangement  as  paneling  P»NN.SYY  where  this  nomenclature  is 
defined  in  Table  VII.  For  example,  paneling  P=45.100 
represents  an  airfoil  with  45  panels  which  are  symmetrically 
arranged  about  the  x  axis  and  which  were  generated  by  equal 
angle  increments  in  the  complex  circle  plane. 

Most  of  the  results  will  be  presented  as  plots  of 
surface  tangential  velocity  error  at  control  points.  This 
procedure  implicitly  assumes  that  the  paneling  is  a  good 
model  of  the  actual  surface.  Referring  to  Figure  43,  the 
exact  and  computed  tangential  velocities  are  determined  at 
the  points  Pg  and  PM  respectively.  The  point  Pg  is 


TABLE  VIII 


Preliminary  Results  for  the  Karman-Tref ftz  Airfoil 


Method 

10 

20 

QC 

21 

30 

31 

32 

Source , 

V 

n 

=  0 

S 

S 

U 

s 

u 

S 

Vortex , 

Vt 

=  0 

s 

S 

u 

s 

u 

u* 

Vortex, 

V 

n 

=  0 

u 

u 

s 

u 

s 

u 

S  =  Successful 
U  =  Unsuccessful 

t  *  =  Non  Oscillatory  Solution, 

But  Incorrect  Lift 

determined  by  requiring  that  it  have  the  same  x  coordinate 
as  PM ' 

Preliminary  Results.  The  method  was  initially  applied 
to  a  19  panel  model  of  the  chosen  Karman-Tref ftz  airfoil. 
Both  source  and  vortex  singularities  were  used  with  the 
singularity  distributions  expanded  about  the  panel  center 
point,  for  various  combinations  of  Q  and  C.  Table  VIII 
indicates  which  of  tt  re  initial  efforts  were  successful. 
Successful  cases  are  those  cases  for  which  at  least  a 
reasonable  solution  was  obtained.  Unsuccessful  cases  are 


those  cases  for  which  the  aerodynamic  influence  matrices 


were  either  algorithmically  singular,  or  for  which  the 
solution  exhibited  an  oscillatory  behavior.  The  asterisked 
case  was  one  for  which  the  solution  seemed  reasonable, 
except  that  the  lift  was  considerably  in  error.  Figures  44 
and  45  show  typical  tangential  velocity  error  results  for 
successful  source  and  vortex  calculations.  For  the  vortex 
case  the  Kutta  condition  was  satisfied  by  placing  an 
internal  point  source  at  Xc  =  . 5  ,  and  specifying  zero 
vorticity  at  the  trailing  edge. 

Several  approaches  were  tried  to  obtain  successful 
solutions  for  all  cases.  For  the  vortex  cases  a  number  of 
alternate  Kutta  conditions  were  used.  These  were  a 
specification  of  zero  vorticity  at  the  trailing  edge  with  an 
internal  point  source,  a  specification  of  net  circulation  in 
an  error  parameter  approach,  and  a  specification  of  zero 
velocity  normal  to  a  trailing  edge  bisector  at  a  point 
slightly  downstream  of  the  trailing  edge  using  both  an 
internal  source,  and  an  error  parameter  approach.  The 
results  of  these  attempts  were  essentially  identical  to 
those  shown  in  Table  VIII.  The  error  parameter  approach  is 
equivalent  to  the  internal  point  source  approach,  based  on 
results  to  date.  The  unsuccessful  cases  were  also  attempted 
using  the  non  centered  form  of  the  series  expansion  with  no 
change  in  the  results. 

For  the  source  case,  Fig  46  shows  an  unsuccessful 


7b.  Effect  of  N  on  Tangential  Velocity  Error,  QC=21,  for  u=.l  vac 
Upper  Surface. 


Figure  47e.  Effect  of  N  on  Tangential  Velocity  Error,  QC=21,  for  u= . 1  rad 
Lower  Surface. 


of  increasing  N  is  to  decrease  the  magnitude  of  the  point 
source  strength  and  thus  the  velocities  induced  by  the 
source  are  correspondingly  reduced.  Figure  48  shows  the 
effect  of  changing  Q  for  a  45  element  case.  This  effect  is 
not  large,  but  it  is  interesting  to  note  that  the  Q=3  cases 
do  not  exhibit  the  source  induced  error  apparent  in  the  Q=2 
cases,  although  the  magnitude  of  the  source  strengths  are  of 
the  same  order.  The  Q=3  case,  however,  has  two  control 
points  on  a  panel,  and  with  this  additional  control  point 
the  normal  flow  on  the  panel  is  more  effectively  controlled. 

Figure  49  shows  the  effect  of  additional  panels  (and 
control  points)  near  the  source  location.  The  airfoil  used 
in  this  figure  is  the  65  element  airfoil  of  Figure  47,  with 
the  addition  of  10  panels  on  each  of  the  upper  and  lower 
surfaces  between  40%  and  60%  chord  for  a  total  paneling  of 
85  elements.  The  point  source  location  remained 
at  Xs=.5  .  The  result  is  that  the  additional  control 
points  in  the  vicinity  of  the  source  control  the  source 
strength  quite  well. 

Effect  of  Point  Source  Location.  A  study  was  conducted 
to  determine  more  precisely  the  effect  of  the  point  source 
location  on  the  solution.  Figures  50a  and  50b  show  the 
tangential  velocity  error  for  a  45  element  airfoil  for  the 
Q*2  and  3  cases  at  a  =  .0  where  the  internal  point  was 
placed  at  one  of  four  locations,  X  .  There  are  two 
interesting  points  to  note  in  these  figures.  The  first  is 
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Figure  41).  Effect  of  N  on  Source  Induced  Tangential  Velocity  Error  for  a  Kannun 
Trefftz  Airfoil. 


Figure  50a.  Effect  of  Point  Source  Location  on  Tangential  Velocity  Error  for 
Karinan-Tref  f  tz  Airfoil,  QC=21. 


that  the  effect  of  changing  from  Q=2  to  Q=3  removes  the 
source  induced  error  for  reasons  discussed  above.  The 
second  is  that  the  source  induced  error  is  indeed  a  very 
local  effect.  Whether  the  source  is  at  X  =.01,  .1  , 

o 

or  .5  ,  the  solution  at  the  trailing  edge  remains 
unaffected,  and  when  the  source  is  at  Xg  =  .9  ,  the 
solution  over  the  leading  80%  of  the  airfoil  is  unchanged. 

The  effects  are  local  because  the  source  strengths  are 
small,  and  their  effect  on  velocity  falls  off  as  the  inverse 
square  of  the  distance  from  the  source.  It  shou-d  be 
mentioned  that  another  way  of  diffusing  the  effect  of  the 
point  source  would  be  to  use  a  distributed  source  on  a  line 
inside  the  airfoil. 

If  the  accuracy  of  the  solution  over  the  whole  airfoil 
is  considered,  these  results  indicate  that  the  best  location 
for  the  point  source  is  very  near  the  leading  edge.  For 
this  case  the  solution  is  excellent  over  99%  o±  the  airfoil, 
while  the  source  induced  error  at  the  nose  is  masked 
somewhat  by  the,  in  general,  larger  error  that  occurs  in 
this  region.  These  errors  are  partly  due  to  the  fact  that  a 
panel  method  will  have  difficulty  accounting  for  the  rapid 
changes  that  occur  in  velocity  as  one  moves  away  from  the 
stagnation  point.  On  physical  grounds  it  might  also  be 
argued  that  for  an  airfoil  with  a  blunt  nose,  a  point  source 
near  the  nose  would  be  able  to  more  effectively  control  the 


oncoming  free  stream  than  would  the  surface  vorticity  in  the 
nose  region. 

Effect  of  Panel  Geometry  Characteristics.  A  study  was 
also  conducted  to  determine  the  effects  of  panel  geometric 
characteristics  on  the  solution.  A  number  of  these 
geometric  parameters  are  involved  in  the  present  method, 
including  the  angle  subtended  by  each  panel,  the  arc  length 
of  each  panel,  and  the  curvature  of  each  panel.  One  would 
intuitively  expect  these  geometric  parameters  to  vary  in  a 
smooth  manner  around  the  airfoil,  making  a  reasonable 
approximation  to  the  actual  airfoil  characteristics.  The 
panel  approximation  is  of  course,  only  a  piecewise 
continuous  representation  of  the  surface,  and  will  exhibit 
discontinuities  of  curvature.  The  curvature  of  a  panel  in  a 
flat  panel  method  is  constant  and  equal  to  zero  over  a 
panel,  while  in  the  present  method  the  curvature  is  constant 
but  not  zero  over  a  panel,  with  variations  from  panel  to 
panel.  One  would  also  expect  the  panel  arc  length  to  vary 
somewhat  smoothly  around  the  airfoil.  One  would  not  expect 
good  results  with  a  very  small  panel  between  two  large 
panels  since  the  small  panel's  control  point  would  be 
overpowered  by  the  nearby  larger  panels.  Several 
investigators  (Henshaw,  Ref  25,  Hess,  Ref  12)  have  suggested 
as  a  rule  of  thumb  that  the  ratio  of  arc  lengths  of  adjacent 
panels  be  less  than  1.5,  but  this  is  merely  a  suggestion 
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which  is  probably  somewhat  dependent  on  the  method  being 
used. 

Figures  51a  -  51c  show  paneling  characteristics  for 

five  different  panel  arrangements  for  the  Karman-Tref ftz 
airfoil  being  considered.  The  angle  subtended  by  each 
panel,  the  panel  radius,  and  a  normalized  panel  arc  length 
are  plotted  against  panel  number  where  panel  one  is  the 
first  panel  at  the  trailing  edge.  Since  the  airfoil  and  all 
these  panelings  are  symmetric  only  upper  surface  quantities 
have  been  plotted.  The  panel  models  compared  to  the  actual 
surface  have  been  plotted  with  an  expanded  vertical  scale  in 
Figures  52a-52d.  As  was  the  case  in  Figure  42  almost  no 
difference  between  the  two  can  be  seen,  except  in  the 
paneling  -  45.101  case  which  will  be  discussed  below.  The 
tangential  velocity  errors  for  these  panelings  are  shown  in 
Figures  53a  and  53b  for  QC=21  and  QC=31  respectively,  with 
a  =  0  radians.  Similar  results  for  ot  =  .  1  radians  are 
shown  in  Figures  53c-53f.  In  these  cases  results  are  shown 
for  both  upper  and  lower  surfaces  since  they  are  not  the 
same  at  non-zero  angle  of  attack. 

The  first  panel  arrangement  which  was  computed  was  the 
p=45.100  scheme,  and  while  the  geometric  characteristics  as 
well  as  the  velocity  error  results  appear  reasonable,  it  was 
felt  that  some  improvements  could  be  obtained,  particularly 
in  the  leading  edge  region.  This  paneling  was  then  modified 
slightly  by  shifting  panels  toward  the  nose  and  the  tail  to 
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Figure  51b.  Paneling  Geometric  Characteristics  for  a  Karman 
Trefftz  Airfoil,  Panel  Radius. 
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Figure  51c.  Paneling  Geometric  Characteristics  for  a  Kartnan 
Trefftz  Airfoil,  Panel  Arc  Length. 
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Figure  52a.  Comparison  of  Panel  Model  with  Actual  Surface 
for  a  Karman-Tref ftz  Airfoil,  Panelingas45 . 101 . 


CHORD  (X/C) 


Figure  52b.  Comparison  of  Panel  Model  with  Actual  Surface 
for  a  Karman-Tref ftz  Airfoil,  Paneling»45 . 102 . 


172 


173 


0.06 


O 

N  ao3 


e 


0.00 

ao  0J2  Ol4  0.6  0.8  LO 

CHORD  (X/C) 

Figure  52d.  Comparison  of  Panel  Model  with  Actual  Surface 
for  a  Karman-Tref ftz  Airfoil,  Panel ing=49. 101. 
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Figure  53d.  Effect  of  Paneling  Character ist ics  on  Tangential  Velocity  Err 
Karman-Tre f f tz  Airfoil,  QC=21.  a^.l  radians.  Lower  Surface. 


VORTEX  DISTRIBUTION  I  I  VORTEX  DISTRIBUTION 

LOWER  SURFACE  I  I  LOWER  SURFACE 


produce  the  P  =  45.101  paneling,  but  in  so  doing,  a  rather 
large  jump  in  all  of  the  panel  geometric  characteristics  is 
introduced  near  the  30%  chord  point  at  panel  No. 14.  The 
effect  on  the  surface  modeling  of  this  modification  can  be 
seen  in  Figure  52a.  Figures  53a-53f  show  the  dramatic 
effect  of  this  discontinuity  on  the  velocity  error,  but  note 
that  the  effect  is  quite  localized  near  the  discontinuity. 

The  P=45.101  paneling  was  then  modified  by  slightly 
moving  slightly  one  of  the  panel  defining  points  for  the 
panel  on  which  the  curvature  jump  takes  place.  The 
resultant  P=45.103  paneling  is  only  slightly  different  than 
the  P=45.101  case,  yet  the  curvature  jump  is  gone  and  the 
resultant  error  plots  are  much  improved  as  shown  in  Figures 
53a-53f . 

The  original  P=45.100  paneling  was  modified  in  a 
different  way  by  removing  panels  at  the  nose  and  tail  while 
requiring  that  the  1.5  rule  hold  in  these  regions,  and 
adding  panels  over  the  mid  section  of  the  airfoil.  This 
arrangement  is  the  P  =  45.102  paneling  and  the  velocity 
error  plots  for  this  case  indeed  show  a  reduction  in  error 
at  the  leading  and  trailing  edges  compared  to  the  P  =  45.100 
case,  indicating  the  validity  of  the  1.5  rule.  This 
paneling  was  then  further  modified  by  adding  panels  at  the 
leading  and  trailing  edge  to  obtain  the  P=49.101  paneling. 
The  errors  for  this  case  are  similar  to  those  for  the 
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P=45.102  arrangement  with  perhaps  slight  improvement  near 
the  leading  edge  which  may  be  attributable  to  the  slightly 
increased  panel  density. 

These  studies  have  shown  that  the  present  method  is 
rather  insensitive  to  changes  in  geometric  paneling 
characteristics  provided  that  these  characteristics  conform 
generally  to  what  might  be  expected  of  a  good  model  of  the 
airfoil;  that  is,  the  characteristics  vary  smoothly  over  the 
airfoil.  The  method  is  sensitive  to  these  characteristics 
only  when  large  changes  occur  over  a  small  part  of  the 
airfoil,  but  in  these  cases  the  sensitivity  effects  are 
localized.  The  effects  are  localized  for  the  same  reason 
that  the  effect  of  the  point  source  was  a  local  one,  i.e. 
the  effect  of  a  panel  decreases  with  distance  from  the 
panel.  The  reason  that  the  method  is  sensitive,  locally,  to 
geometric  characteristics  is  that  the  method  is  capable  of 
modeling  geometry  very  well,  particularly  as  N  is  increased. 
This  capability  for  faithful  geometric  representation  is  one 
of  the  attractive  features  of  panel  methods  in  general. 

Effect  of  Panel  Curvature  A  study  was  also  conducted 
to  determine  the  effects  of  panel  curvature  on  the  computed 
solution.  The  parameter  controls  the  curvature  of  each 
panel  by  moving  the  third  point,  through  which  the  circle  is 
drawn,  in  from  the  originally  specified  point  on  the  airfoil 
surface  toward  the  straight  line  connecting  the  panel  end 
points.  Thus  8=0  would  be,  in  the  limit,  a  flat  panel, 
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and  8=1  is  the  normal  case  where  the  circle  is  drawn 
through  the  three  input  points  on  the  airfoil  surface. 

Figures  54a-54f  show  results  for  the  P=45.1Q2  airfoil 
with  a  point  source  at  Xc  =  .01  for  a  range  of  B's  from 
1  to  0.2,  and  for  two  angles  of  attack.  Note  that  the 
errors  increase  as  the  panel  becomes  flatter.  Also 
note  that  as  the  panel  becomes  flatter  the  error  is  larger 
for  the  three  term  series  than  for  the  two  term  series  in 
all  cases.  It  seems  that  as  the  panel  becomes  flatter,  or 
approaches  a  lower  order  representation  of  the  surface,  a 
lower  order  representation  of  the  singularity  (i.e.  a  2  term 
rather  than  a  3  term  series)  is  sufficient  to  produce  a 
certain  level  of  error.  Of  course  these  errors  are  larger 
than  are  obtained  by  using  the  full  circular  arc  panel.  The 
reason  that  better  results  are  obtained  using  the  more 
highly  curved  panels  is  that  they  provide  a  better 
representation  of  the  surface  in  terms  of  location  of  the 
singularity,  and  location  of  the  control  point. 

Effect  of  Angle  of  Attack.  The  results  presented  to 
this  point  have  included  angles  of  attack  ot=.0  ,  and  .1 
radians,  and  have  exhibited  no  strong  sensitivity  to  the 
angle  of  attack.  To  study  more  carefully  the  effects  of 
angle  of  attack,  an  N=19  element  symmetrically  paneled 
airfoil  was  computed  at  various  a's  between  0.,  and  .55 
radians  for  both  the  two  and  the  three  term  singularity 
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Figure  54f.  Effect  of  Panel  Curvature  on  Tangential  Velocity  Error  for 
Trefftz  Airfoil,  QC=31,  u=.l  radians,  Lower  Surface. 


series.  The  results  of  these  calculations  are  presented  in 

Figure  55  which  shows  excellent  agreement  for  the  lift 

coefficient  over  the  entire  angle  of  attack  range  for  both 

two  and  three  term  cases.  Figure  56  shows  the  error  in 

lift  coefficient;  that  is  CL  “  CL  “  CL  .  Note 

er  ex  comp 

that  both  the  two  and  three  term  curves  are  linear  with  a  , 

and  that  the  error  is  quite  small.  Also,  the  error  for  the 

three  term  case  is  less  than  for  the  two  term  case,  and  the 

slope  of  the  three  term  curve  is  less  than  that  for  the  two 

term  case  so  that  the  difference  increases  with  a  .  A 

consequence  of  the  linearity  of  these  curves  is  that  the  fig 
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relative  error,  defined  as  L  ,  comP  .  is 

rel  c 
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essentially  a  constant.  For  _  _  _  2  ,  and 
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for  Q=3 ,  CLrel=  .00578 


Q=2,  CLrel=. 00697 


Effect  of  Control  Point  Location. 


important 


question  which  arises  in  the  discussion  of  any  panel  method 
is  the  sensitivity  of  the  solution  to  control  point 
placement.  In  the  present  method,  the  control  point 
location  is  assummed  to  be  on  the  circular  arc  which 
represents  the  surface,  and  similarly,  the  normal  to  the 
surface  is  represented  by  the  normal  to  the  circular  arc. 
Thus  there  is  an  error  in  both  the  control  point  location, 
and  the  surface  slope.  Since  surface  location  and  slope  can 
be  computed  exactly  for  a  Karman-Tref f tz  airfoil,  a  study 
was  done  to  determine  the  sensitivity  of  the  current  method 


to  these  errors.  Since  these  two  parameters  are  independent 


190 


four  cases  can  be  considered: 

1.  the  location  and  slope  are  computed,  as  in  the 
basic  method,  using  circular  arc  panels. 

2.  the  location  is  exact  but  the  slope  is  computed 
from  the  circular  arc. 

3.  the  location  is  computed  from  the  circular  arc 
panels  but  the  slope  is  exact,  and 

4.  both  the  location  and  the  slope  are  exact. 

For  the  two  term  model,  in  which  there  is  one  control  point 
per  panel,  the  third  panel  defining  point  is  used  as  the 
control  point  in  cases  two  and  four.  For  the  three  term 
model,  where  two  control  points  are  required  on  a  panel 
instead  of  one,  this  approach  cannot  be  used;  so  for  the 
three  term  model  only  cases  one  and  three  have  been  studied. 

Figures  57a-e  and  58a-e  show  tangential  velocity  error 
results  for  the  four  cases  as  a  function  of  several 
parameters.  Figures  57a  and  57b  are  results  for  a  nine 
element  airfoil  at  a  =  0  for  QC=21  and  QC=31  .  While 

these  figures  do  not  give  a  clear  indication  of  which 
combinations  of  slope  and  control  point  parameters  are 
superior,  they  do  show  that  contrary  to  what  might  be 
expected,  the  use  of  the  exact  slope  and  control  point 
location  is  clearly  not  the  best.  Results  for  a  *  . 1 
radians  shown  in  Figures  57c-e  are  similar.  However,  results 
for  a  ten  element  airfoil,  given  in  Figs  58a-e,  show  that 
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Figure  57a.  Effect  of  Control  Point  Location  and  Slope  on 
Tangential  Velocity  Error  for  a  Karman-Tref ftz 
Airfoil,  N=*9,  QCa2l,  a=.0  radians. 


AD-A124  896 
UNCLASSIFIED 


COMPUTATION  OF  INCOMPRESSIBLE  POTENTIAL  FLOU  OVER  AN 
AIRFOIL  USING  A' HIGH.  .  <U>  AIR  FORCE  INST  OF  TECH 
URIGHT-PATTERSON  AFB  OH  SCHOOL  OF  ENG I.  .  J  DEJONGH 
HUG  82  AFIT/DS/HA/82-i  F/G  12/1 


CHORD  (X/C) 


Figure  57b.  Effect  of  Control  Point  Location  and  Slope  on 
Tangential  Velocity  Error  for  a  Karraan-Tref ftz 
Airfoil,  N“9,  QC=*31,  a=».0  radians. 
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Figure  57c. 
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Effect  of  Control  Point  Location  and  Slope  on 
Tangential  Velocity  Error  for  a  Karman-Tref ftz 
Airfoil,  N=9,  QC=21,  a=.l  radians,  Upper  Surface 


TANGENTIAL  VELOCITY  ERROR 


\  \  / 

\  x  . . 

\  / 


a4 

CHORD  (X/C) 


Figure  57d.  Effect  of  Control  Point  Location  and  Slope  on 
Tangential  Velocity  Error  for  a  Karman-Tref ftz 
Airfoil,  QC-21,  a= .  1  radians,  Lower  Surface 
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Figure  58c.  Effect  of  Control  Point  Location  and  Slope  on 
Tangential  Velocity  Error  for  a  Karman-Tref f tz 
Airfoil,  N*10, QC»21, a= . lradians, Upper  Surface. 
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Figure  58d.  Effect  of  Control  Point  Location  and  Slope  on 
Tangential  Velocity  Error  for  a  Karman-Tref ftz 
Airfoil,  N“10 , QC*21 , a* . 1  radians , Lower  Surface 
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Figure  58e.  Effect  of  Control  Point  Location  and  Slope  on 
Tangential  Velocity  Error  for  a  Karman-Tref ftz 
Airfoil,  N»10,  QC-31,  ct«.l  radians. 
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the  error  using  the  basic  method  is  generally  small  and 
relatively  constant  over  the  center  part  of  the  airfoil, 
while  the  other  cases  give  somewhat  larger  errors  which  also 
tend  to  vary  more  dramatically.  As  before,  the  use  of  the 
exact  slope  and  control  point  location  clearly  does  not  lead 
to  the  best  solution. 

It  is  instructive  to  look  at  the  slope  percentage  error 
that  results  from  using  a  computed  control  point.  Figure  59 
shows  the  slope  error  for  a  19  element  and  a  45  element 
airfoil  with  QC=21  and  Q031,  using  a  computed  control  point 
(i.e.  the  basic  method).  Note  that  the  slope  error  is 
constant  and  small  over  most  of  the  airfoil  for  both  element 
numbers  and  for  both  values  of  Q.  It  is  also  clear  that 
increased  panel  density  in  the  nose  region  significantly 
reduces  the  slope  error  there.  Since  the  errors  on  one  or 
two  panels  near  the  leading  edge  are  relatively  large 
compared  to  the  rest  of  the  airfoil,  one  might  suppose  that 
these  errors  account  in  part  for  the  relatively  larger 
errors  in  tangential  velocity  that  have  been  noted  in  the 
nose  region.  This  could  not  be  the  only  cause  of  these 
errors,  though,  since  larger  errors  are  encountered  in  the 
trailing  edge  region  as  well,  yet  the  slope  errors  in  this 
region  are  very  small. 

This  study  has  shown  that,  contrary  to  what  might  be 
expected,  use  of  exact  control  point  location  and  slope 
information  in  the  present  method  does  not  lead  to  improved 


accuracy  in  the  solution.  Since  the  representation  of  the 
surface  geometry  in  terms  of  slope  error  was  seen  to  be  very 
good,  other  sources  of  error  in  the  method  probably  drive 
the  error  in  the  solution.  The  most  likely  sources  of  error 
are  the  discrete  application  of  the  boundary  conditions,  and 
the  series  approximation  to  the  singularity  distribution. 
Thus  the  additional  input  data  which  would  be  required  to 
use  exact  slope  and  control  point  location  is  not  justified, 
and  the  use  of  computed  slope  and  location  information  is 
completely  adequate  to  represent  the  surface  geometry,  at 
least  as  far  as  the  present  method  is  concerned. 

Effect  of  Camber.  The  airfoils  which  have  been  studied 
up  to  this  point  have  been  symmetric.  Now,  the  present 
method  will  be  applied  to  a  slightly  cambered  Karman-Tref ftz 
airfoil  which  has  a  camber  parameter  of  3ir/4  ,  a  trailing 
edge  angle  of  0.2356  radians,  and  a  zero  lift  angle  of 
attack  of  -.0275  radians.  The  airfoil  was  paneled  with  45 
elements  equally  spaced  in  the  circle  plane,  which  produces 
a  slightly  asymmetric  paneling  (designated  as  P=45.201)  in 
the  airfoil  plane.  The  basic  method  was  used 
with  Xs=.01  and  8=1  .  Results  for  the  QC=21  case  at 
angles  of  attack  0.0  and  0.1  radians  are  shown  in  Figures 
60a  and  60b.  Both  upper  and  lower  surfaces  are  plotted 
since  the  paneling  is  not  symmetric.  Figures  60c  and  60d 
show  similar  results  for  the  QC=31  cases.  These  curves 
exhibit  characteristics  similar  to  those  noted  earlier  for 
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Tungential  Velocity  Error  for  a  Cambered  Karman-Tre  f  ftss  Airfoil 
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the  symmetric  airfoil.  These  velocity  errors  near  the  nose 
do  appear  to  be  larger,  but  this  may  be  due  in  part  to  the 
fact  that  the  point  source  was  placed  very  near  the  nose, 
and  also  to  the  fact  that  the  stagnation  point  at  the  nose 
is  not  near  the  control  point  of  the  nose  panel.  For  the 
symmetric  airfoil  cases  with  a  similar  paneling  the  nose 
panel  control  point  was  located  at  the  stagnation  point,  so 
that  the  correct  solution  was  obtained  at  that  point. 

Summary 

In  this  chapter  the  method  of  circular  arc  panels  was 
applied  to  several  different  types  of  airfoils,  and  the 
characteristics  which  define  the  method  were  varied 
systematically  to  determine  their  effect  on  solution 
accuracy.  The  method  was  first  applied  to  a  Joukowski 
airfoil,  an  NACA  0024  airfoil, a  thin  symmetric  airfoil,  and 
a  Karman-Tref ftz  airfoil  using  different  combinations  of 
singularity  type  and  varying  the  number  of  terms  in  the 
series  and  the  degree  of  continuity  imposed.  These 
preliminary  studies  showed  that  accurate  results  were 
consistently  obtained  for  different  types  of  airfoils  and 
for  lifting  and  non-lifting  cases  by  using  the  following 
approach  (which  can  be  compared  with  the  characteristics 
shown  in  Table  I):  a  2  or  3  term  (linear  or  quadratic) 
vortex  distribution  is  p  laced  on  each  panel,  and  continuity 
of  the  distribution  is  enforced  at  panel  junctures.  The 
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panels  are  piecewise  continuous  circular  arc  elements 
generated  from  the  surface  geometry  with  no  series  expansion 
approximations.  The  boundary  condition  imposed  at  control 
point  is  zero  external  normal  velocity.  The  Kutta  condition 
is  met  by  specifying  zero  vorticity  at  the  trailing  edge  on 
both  upper  and  lower  surfaces,  and  an  internal  point  source 
is  added  to  close  the  formulation.  All  integrations  are 
performed  analytically  for  maximum  computational  efficiency. 
This  basic  method  was  then  exercised  on  both  symmetric  and 
cambered  Karman-Tref ftz  airfoils  at  different  angles  fo 
attack  to  determine  the  effects  of  N,  Q,  panel  geometry, 
point  source  location,  panel  curvature,  and  controfl  point 
characteristics  on  the  accuracy  of  the  method.  The 
following  conclusions  can  be  drawn: 

1.  Increasing  N  and/or  Q  produces  more  accurate 
results. 

2.  The  method  is  somewhat  sensitive  to  paneling 
geometric  characteristics  (panel  subtended  angle,  panel 
radius,  and  panel  arc  length),  but  the  effect  is  local. 

3.  The  effect  of  the  internal  point  source  can  be 
relatively  large,  but  it  is  very  localized  and  can  be 
controlled  by  using  additional  control  points  in  the 
vicinity  of  the  source. 

4.  The  accuracy  of  the  solution  generally  increases 
as  the  curvature  of  the  panels  is  varied  from  nearly  flat  to 
the  circular  arc  model. 


5.  The  effect  of  exact  representation  of  control 
point  location  and  slope  does  not  lead  to  more  accurate 
solutions  compared  to  results  based  on  the  computation  of 
the  location  and  slope  from  the  circular  arc  model. 

6.  The  method  produces  good  results  over  a  range  of 
angles  of  attack,  although  the  accuracy  decreases  linearly 
as  angle  of  attack  increases. 

The  next  chapter  will  summarize  the  development  of  the 
present  method,  draw  conclusions  concerning  the  application 
of  the  method  to  the  circular  cylinder  and  to  airfoils,  and 
present  suggestions  for  future  work  in  this  area. 


VI .  Conclusions  and  Recommendations 

Conclusions 

The  purpose  of  this  effort  was  to  develop  a  better 
understanding  of  the  effects  of  the  several  characteristics 
involved  in  a  panel  method  solution,  and  to  provide  guidance 
and  understanding  for  the  further  development  of  two  and 
three  dimensional  panel  methods.  To  reach  this  goal,  a  new 
panel  method,  based  on  the  fundamental  concepts  of  potential 
theory  and  on  a  simple  approach  to  curve  approximation,  has 
been  developed.  This  method  used  a  new  approximating 
element,  the  circular  arc;  and  a  new  singularity 
representation,  the  sine  series. 

The  method  was  initially  applied  to  the  problem  of  flow 
over  a  circular  cylinder  and  the  effects  of  varying  several 
parameters  were  studied.  This  effort  showed  that  the 
current  method  was  capable  of  accurate  results  (which  were 
noted  earlier),  and  it  allowed  an  assessment  of  the  effects 
of  the  many  characteristics  which  impact  a  solution.  This 
assessment  was  used  to  develop  the  method  further. 

Based  on  these  studies  of  the  circular  cylinder  the 
method  was  applied  to  a  Joukowski  airfoil,  an  NACA  0024 
airfoil,  a  thin  symmetric  airfoil  and  a  Karman-Tref ftz 
airfoil.  Initial  studies  were  performed  to  assess  the 


”■* ’  applicability  of  the  method  to  airfoil  shapes  as  a  function 

of  singularity  type,  number  of  terms  in  the  series  (Q), 
degree  of  continuity  (C),  type  of  boundary  condition,  and 
Kutta  condition  formulation.  Source  and  vortex 

distributions  were  used,  while  values  of  Q  varied  from  1  to 
3  and  values  of  C  varied  from  0  to  2.  For  the  vortex  cases 
3  types  of  Kutta  condition  were  investigated:  an  error 

parameter  approach  with  a  trailing  edge  bisector  condition, 
an  internal  point  source  with  a  trailing  edge  bisector 
condition,  and  an  internal  point  source  with  a  specification 
of  zero  vorticity  at  the  trailing  edge.  These  studies 
indicated  that  the  method  was  not  sensitive  to  the  type  of 
jP*  Kutta  condition  used.  Also  not  all  combinations  of  Q  and  C 

yielded  acceptable  solutions,  depending  on  the  type  of 
singularity  which  was  used. 

As  a  result  of  these  preliminary  studies  a  basic  method 
was  chosen  for  further  investigation.  This  method  used  a 
continuous  2  or  3  term  vorticity  distribution  with  a  normal 
velocity  boundary  condition,  and  an  internal  point  source 
with  zero  vorticity  at  the  trailing  edge  to  satisfy  the 
Kutta  condition.  This  basic  method  was  then  applied  to  a 
symmetric  Karman-Tref ftz  airfoil  and  detailed  studies  were 
conducted  to  determine  the  effects  of  number  of  panels  (N), 
number  of  terms,  point  source  location,  geometric  paneling 
characteristics,  panel  curvature,  angle  of  attack,  control 
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point  location  and  slope,  and  airfoil  camber.  Conclusions 
which  can  be  drawn  from  this  study  are: 

1.  The  method  produces  very  accurate  solutions  over 
the  major  part  of  the  airfoil,  with  the  largest  errors 
occuring  at  the  leading  and  trailing  edges  (i.e.  the 
stagnation  point)  .  These  errors  are,  however,  always  small 
compared  to  the  free  stream  and  are  small  compared  to  the 
exact  solution  except  at  points  next  to  the  stagnation 
points. 

2.  Significant  error  reduction  occurs  as  N  is 
increased,  and  reasonable,  though  not  as  large  a  reduction, 
occurs  as  Q  is  increased  from  2  to  3. 

3.  The  effect  of  point  source  location  is  large  but 
is  very  local.  It  was  found  that  these  source  induced 
errors  can  be  effectively  controlled  by  either  a  3  term 
series,  or  by  placing  additional  control  points  near  the 
source. 

4.  The  method  is  generally  insensitive  to  minor 
variations  in  paneling  as  long  as  the  geometric  parameters 
governing  the  paneling:  That  is,  panel  curvature  and  panel 
subtended  angle,  vary  in  a  smooth  manner  around  the  airfoil. 
Additionally,  the  requirement  that  adjacent  panels  maintain 
a  3:2  or  less  ratio  in  arc  length  was  found  to  be  effective 
in  reducing  errors,  particularly  at  the  trailing  edge. 

5.  The  effect  of  panel  curvature  is  that  accuracy 
increases  as  the  curvature  increases  from  that  of  a  nearly 


flat  panel  to  that  of  a  circular  arc  panel  with  three  points 
on  the  circular  arc  coincident  with  the  airfoil  surface  for 
both  the  2  and  3  term  series  expansions. 

6.  The  accuracy  of  the  solution  decreases  slightly  as 
the  angle  of  attack  increases. 

7.  It  was  found  that  improved  accuracy  is  not 
generally  obtained  when  either  exact  control  point  location 
or  slope  information  is  used  as  opposed  to  when  these 
quantities  are  computed  from  the  circular  arc  panel.  Since 
one  might  expect  that  the  additional  exact  information  would 
improve  the  solution,  the  fact  that  it  does  not  indicates 
that  other  errors  inherent  in  the  formulation,  such  as  the 
singularity  formulation  and  the  discretization  process 
itself,  may  be  the  primary  causes  of  error. 

8.  The  method  provides  accurate  results  for  a  non- 
symmetric  airfoil,  although  the  accuracy  is  not  quite  as 
good  as  for  the  symmetric  case. 

Recommendat ions 

Several  areas  for  further  work  have  become  apparent 
during  the  course  of  this  research.  These  include  improving 
the  method  as  a  two  dimensional  tool  and  extending  the 
method  to  the  three  dimensional  case.  In  terms  of 
improvement  of  the  two  dimensional  method,  further  study  to 
improve  the  solution  at  the  leading  and  trailing  edges  could 
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be  undertaken.  A  possible  approach  to  doing  this  would  be 
to  modify  the  singularity  distribution  on  one  or  several 
panels  at  the  leading  or  trailing  edges.  The  rationale  for 
this  approach  is  that  the  gradients  in  singularity  strength 
are  largest  in  these  regions,  and  it  is  possible  that  the 
same  series  representation  of  the  singularity  can  not 
adequately  model  these  gradients.  For  example,  one  might 
use,  at  the  trailing  edges,  a  singularity  strength  which  is 
proportional  to  the  square  root  of  the  arc  length  measured 
from  the  trailing  edge  because  this  singularity  will  go  to 
zero  at  the  trailing  edge  more  quickly  than  will  a  linear 
function.  A  disadvantage,  however,  to  using  more 
complicated  representations  of  the  singularity  is  that 
numerical  integration  might  be  required  to  obtain  the 
influence  coefficients  for  the  panels  involved. 

Another  approach  to  quickly  improving  solution  accuracy 
is  based  on  the  observation  that  the  tangential  velocity 
errors  near  the  leading  and  trailing  edges  do  not  vary 
smoothly.  This  is  true  over  a  larger  portion  of  the  airfoil 
for  the  3  term  series  expansion  cases  as  well.  As  a  way  of 
smoothing  these  curves  and  reducing  the  overall  error  in  the 
solution  a  new  velocity  curve  could  be  fitted  to  the 
calculated  results  using  a  least  squares  procedure,  or  some 
type  of  averaging  procedure.  A  systematic  study  of  a 
particular  algorithm  for  doing  this  would  be  required  to 
establish  the  validity  of  this  approach  for  the  case  of  a 


general  airfoil. 

A  third  area  for  additional  study  is  related  to  the 
initial  results  presented  in  Chapter  V.  Further 
investigation  of  the  failure  of  certain  parameter 
combinations  to  yield  solutions  is  required  to  fully 
understand  the  proper  way  of  to  numerically  solve  the 
integral  equations  of  potential  flow  using  the  panel  method 
approach.  This  should  include  further  study  of  the 
application  of  the  Kutta  condition,  particularly  in  the  case 
in  which  a  solution  was  obtained  using  the  internal  velocity 
boundary  condition. 

The  method  could  also  be  extended  to  three  dimensions. 
For  example,  consider  the  case  of  a  finite  wing.  Paneling 
in  the  spanwise  direction  would  have  to  be  developed.  A 
scheme  using  flat  panels  spanwise,  or  one  using  curved 
panels  whose  radii  varied  in  the  spanwise  direction  could  be 
investigated.  Another  alternative  would  be  a  complete 
numerical  integration  in  the  spanwise  direction  coupled  with 
analytic  intergration  chordwise.  If  circular  arc  panels 
were  developed  spanwise,  the  resulting  integrations  to 
obtain  the  influence  coefficients  would  probably  have  to  be 
performed  numerically  since  the  resulting  equations  are 
elliptic.  Another  approach  would  be  to  use  circular  arc 
panels  only  at  the  leading  edge  or  at  the  wing  tips  since 
these  are  the  regions  where  wing  surfaces  typically  exhibit 
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the  largest  curvatures. 

An  additional  extension  of  the  two  dimensional  method 
which  would  have  application  in  three  dimensional  cases  as 
well  would  be  to  develop  a  procedure  for  computing  forces 
and  moments  on  the  airfoil  or  wing  using  computed  velocities 
and  the  surface  paneling.  This  is  an  area  where  the  details 
of  the  surface  model  could  play  an  important  role,  and  an 
assessment  of  the  effect  of  the  circular  arc  panel  model  on 
these  quantities  would  be  a  worthwhile  result. 
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Appendix  A 


Continuity  Coefficient  Matrices 


The  matrix  equations  governing  singularity  strength 
continuity  across  panel  junctures  are  given  by  eq.51  as 


EMkI  -  °  c - o-1'---3 


(51) 


where  the  C^  are  obtained  from  eqs.49  and  50.  These 
matrices  have  the  same  form,  shown  in  eq.70,  for  all  values 
of  C  and  k.  Note  that  the  elements  of  the  matrices  do  de¬ 
pend  on  C  and  k  but  they  have  not  been  marked  as  such  to  re 
duce  the  complexity  of  the  notation.  For  a  given  C  and  k 
these  elements  depend  on  panel  geometric  characteristics  as 
defined  in  Fig  12.  The  elements  themselves  are  defined  in 
Tables  IX  and  X  for  C=l,2  and  3.  For  the  case  C=0  all  el¬ 
ements  are  zero.  It  should  be  noted  also  that  the  last  row 
in  each  matrix,  which  destroys  its  bandedness,  defines  a  re 
lationship  between  the  first  and  the  last  panel.  In  some 
cases,  such  as  where  a  slope  discontinuity  exists  across 
these  panels,  the  conditions  of  continuity  may  be  slightly 


modified. 


Appendix  B 


Velocity  Influence  Coefficients 

In  this  appendix  the  normal  and  tangential  velocity 
influence  coefficient  matrices  [Rfe]  and  [Tk]  wiH  be 
computed.  These  matrices  give  the  velocity  at  any  point 
in  the  field  that  is  induced  by  a  source  or  vortex  singu¬ 
larity  distributed  on  a  circular  arc.  Also  the  velocity 
influence  coefficients  for  a  reverse  curvature  panel,  which 
are  needed  to  model  a  general  airfoil,  will  be  developed. 

Source  Distribution 

Referring  to  Fig  12  and  61,  the  problem  is  to  compute 
the  velocity  at  a  point  P(r,0)  due  to  a  source  distribu¬ 
tion  on  a  circular  arc  panel.  The  velocity  is  given  by 

?<r,e)  -5^*  ■  sll  V  ?  if  ;s)  <71> 

where 

<J>  53  J a  log  R  dl 
B 

R2  =  r2+a2-2ar  cos(9-0o) 

Now 


B 
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where 


Thus 


%  R 

dl  =  ad8  o  ,  -jrp  =  r  -  a  cos(8-0o) 


or 


and 


M+6 


34> 

3r 


f  0(0.)  ■  .T-a  eos(e-e,) - 

r2+a2-2ar  cos(8-e0) 


1-6 


Now  let 


Then 


d  =  L  D  =  r2+a2 
a  ’  2ar 


3(j) 

3r 


=  a 


V6 

/ 

9m's 


g(9n)[r-a  cos(8-8 o )]d8 q 
.  2arLD-cos( 8-0  0 )  J 


34) 

3r 
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2ar 


0M+6 
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a( 8  o ) [d-cos( 9-8  q ) ]d8  o 
D— cos ( 9—8 o )  ” 


VS 


V  (r  0) 

J*'  )  ^ 


=  -L-  /* 

4ird  J 
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q(8  o ) [d-cos( 8-8  0 )]d8 Q 
D-cos(0-0o) 


Similarly  for  : 


where 


3R  _  ar  sin(9-90) 
39  R 


Therefore 


ve(r,9) 


1 

4itd 


eM+6 


a(9 o )sin( 


D-cos( 9 


Now 


Q-l 

cr(  0  o  )  -  ^  qk  sink(90-9M) 
k=0 


so  that  the  general  term  is 


9-9  o )d9  o 

■=07) 


Integration  for  rfa  .  To  integrate  these  expressions 
the  following  transformation  is  made: 

Let 


x  =  0n-e 


Xi  =  V<$-e  ,  X2  =  9M+6-e 


M 


9  o  —9,,  =  X+9-0 


M 


M 


Therefore 


Vrk(rf0) 


A  2  i_ 

qk  y  sin  (X+0-0u)fd-cosXldX 


4ird 


D-cosX 


X, 


Vflk(r,0) 


k  f 
irFcT  J 

Xx 


x2  _J_ki 


-q,,  sin  (X+0-0M)sinXdX 


D-cosX 


(74) 


(75) 


Now  let 

W  *  D-cos  X 

and 

y  =  l-cos2X  =  l-(D-W) 2  =  ( 1-D2 )+2DW-W2 
or 

y  ■  a+bw+cw2 

where 

a  =  1-D2 
b  =  2D 
c  *  -1 
Note  that 


dw  =  sinX  dX 


The  integrals  in  eqs.74  and  75  can  be  obtained  by  reduction 
to  the  following  standard  forms  (Refs  69  and  70): 


sin” 1 ( cosX) 


-sin” 1 (D-W) 


2 

/D2-l 


Tan 


I  D-l 


+  ( 1-D2 


f  ■  -/y  +  D  /*— 

J  /y  ^  /y 


f  widw  „  ^W-SD  ^  +  (D2+i)  rjL 

J  fy  *  /y  w 


Note  that  eq.76a  evaluated  at  the  integration  limits 
Xz 

I  2Z  »  -[sin“1(cosX2)-sin“l(cosXl)]  *  X2-X 


(76a) 


(76b) 


(76c) 


(76d) 


(76e) 


gives 


Using  the  above  in  eqs.74  and  75,  and  letting 

Tan 


m  _-M  /D*7-!  m...  eM-6"9 
-  Tan  j  D_1  Tan  2 


(77) 


F2 (r  ^9 )  =  In 


D-cos( 5-0+0m)  / 

D-cos(-6-9+9m)  ( 


(78) 


X2=6-0+9m  Xl— 6-0+0m  (o-6+0-0m 

The  following  expressions  for  the  velocities  in  the  r  and 
0  directions  due  to  a  source  distribution  on  a  circular 
arc  are  obtained  for  k=0  to  k*=3 
For  k=0 

Vr 0 (r  ,0 )  =  {26+(d-D)Fi (r  ,0 )}  (79a) 

v9°(r,0)  *  Sd  Fz(r >9)  (79b) 

For  k=l 

V^r,9)  -  ^{[(d-D)F2(r  ,0)-cos(X2)+cos(X1  )]cos(0-eM) 


+[sin(X2 )-sin(Xj )-(d-D)2S+DF! (r  ,9 ) ]sin( 0-0M) } 


V(r>e)  “  i^^2D5+sin(x2)-sin(X1)-(D2-l)FJ(r>0)]cos(9-eM) 


+[DF2(r>6)+cos(X2 )-cos(X1 ) ]sin( 8-0M) }  (80b) 
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For  k=2 
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and  for  k=3 
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*(sinX2-sinX i )-  ^(sin2X2-sin2X ! ) 

+  6D( 3-2D2  )+  2  F1(r,0)j+3sin(0-eM)cos2(0-eM) 

■  j^D(l-D2)F2(r,  0)-(3D2-l)(cosX1-cosX1) 

-  ~  [(D-cosX2  )  2-(D-cosX1  ) 2  ]+  |[(D-cosX2)3 
-(D-cosXi ) 3jj  +3sin2( 0-0M)cos( 0-0^) [D2(sinX2-sinXj ) 
+26D(D2-i)+  ^  [(D-cosX2 )sinX2-(D-cosXi )sinXj ] 

-  sln»X,-sln»Xi  +  Dj(l-p  F,  (r,  9 )  ]+sin  >  ( e-eH) 


•  |^D2F2  (r.  0  )+3D2  (cosX2-cosXj )+  jH[(D-cosX2  )2 
-(D-cosX! )2]-  i[ (D-cosX2 ) 3-(D-cosXi ) 3 ]] }  (83b) 


Limit  as  r-»-a  .  The  above  formulas  represent  the  ef¬ 

fect  of  a  source  distribution  on  a  circular  arc  at  a  point 
PCr^e)  .  They  are  clearly  valid  for  r^a  ,  but  if  P  lies 
on  the  arc  itself  they  exhibit  special  limiting  behavior. 

In  this  case  P(r  0)-*-P(a  0)  with  0e  ( 0,.-<$ ,  0,,+S )  and  the 
unique  limiting  behavior  is  contained  in  terms  with  the  fac¬ 
tor  Fi(r,0)  .  Now  consider  the  factor  /ti2-i/(D-l)  in 

Eq.77.  As 

r^a  '  Torn  -  00  ' 

But  the  tangent  term  multiplying  this  factor  changes  sign  de¬ 
pending  on  the  value  of  0  relative  to  the  panel.  Refer¬ 
ring  to  Fig.  62  consider  a  point  Pi  such  that  0 1£(0  _<S  f  0„+6 ) . 
Then  both  tangent  terms  in  Eq.77  are  negative;  i.e. 

0„+6-0i  Q^.-6-Qi 

Tan  -  <0  and  Tan  -2-g -  <0 

Thus  as  r+a  with  r>a  the  inverse  tangent  difference 
term  in  Eq.77  becomes 

Tan” 1  ( -°° )  -  Tan- 1  ( -®)  *  0 

The  same  is  true  for  the  case  of  P3  ,  except  that  the 
tangents  are  positive.  For  the  case  of  P2  ,  however, 


and 


0m-6-62 


0 


Thus 


Tan”1  (+«>)-Tan-1  (-<»)  =  J  -(-  J)  =  tt 


0 


As  an  example,  consider  Vr(r,9)  for  the  k=0  case,  and 
let  r-*-a+  .  Note  that 


lim 

r-*-a 


d-D 


sgn  (r-a) 


Therefore 


+ 

as  r+a 


d-D 


1 


Thus  for  the  case  corresponding  to  P2  in  Fig  61 

vr(a+»  9  )  *  [26+2tt ] 

While  for  the  cases  corresponding  to  Pj  or  P3  , 

Vr(a+,0)  -  3£[26] 

An  interesting  feature  of  this  expression  is  that  for  the 
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1 


k=0  case  the  dependence  of  9  is  removed  by  the  limiting 
process. 

Transformation  to  Global  Coordinates.  Eqs  79-83  rep¬ 
resent  velocities  at  a  point  P  in  a  polar  coordinate  sys¬ 
tem  based  on  the  circular  arc  panel.  To  obtain  the  global 
influence  coefficient  matrices  [  ]  and  [  T^  ]  these 

equations  must  be  transformed  by  a  simple  rotation  into 
tangent  and  normal  velocities  defined  at  the  point  P 

Reverse  Curvature  Panels 

The  formulas  that  have  been  developed  are  written  in 
a  polar  coordinate  system  in  which  the  equation  of  the 
panel  is  r=>a  .  In  addition,  the  direction  of  increasing 
arc  length  along  an  element  is  assumed  to  be  counterclock¬ 
wise  in  such  a  coordinate  system.  This  type  of  panel  will 
be  referred  to  as  a  standard  panel.  A  general  cambered 
airfoil,  however,  will  have  regions  of  curvature  opposite 
to  that  of  a  standard  panel.  By  applying  a  transformation 
the  formulas  for  a  standard  panel  can  be  used  to  compute 
the  velocities  due  to  a  reverse  curvature  panel. 

The  problem  is  to  compute  the  velocity  at  a  point  P 
due  to  a  distribution  on  arc  B  (Fig  63)..  Let  $  be  a 
vector  defined  by  connecting  the  panel  endpoints  Pj  and 
P2  ;  and  let  fi  be  the  direction  of  the  normal  to  T 

A  A 

The  point  P  is  then  obtained  by  reflecting  P  about  T. 

A  A 

Now  it  is  clear  that  the  velocity  at  P  in  the  T 
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direction  due  to  arc  A  equals  the  velocity  at  P  in  the 


T  direction  due  to  arc  B  .  Similarly  the  velocity  at  P 
in  the  fi  direction  due  to  arc  A  can  be  related  to  the 


equivalent  velocity  at  P  due  to  arc  B  .  That  is, 


V-  (P)  =  V-  (P) 
XB  1 A 


Vn  (P>  = 

nB  nA 


Thus  the  velocity  at  a  point  P  due  to  either  a  standard 
or  a  reverse  curvature  panel  can  be  computed  in  terms  of 
the  local  panel  coordinate  system.  The  components  are  then 
transformed  into  a  global  system  in  which  the  point  is 
specified. 


Vorticity  Distribution 

In  two  dimensions  the  velocity  components  induced  by  a 

vortex  distribution  are  directly  related  to  those  developed 

above  for  a  source  distribution.  If  the  velocity  due  to  a 

source  is  ^  where 
s 


rs  ‘  Vr  er+V0  < 
s  s 


Then  the  velocity  due  to  a  vortex  f  is 


^  =  -V  e  e. 

v  es  r  rg  0 

Thus  the  potential  flow  problem  can  be  immediately  solved 
for  a  vorticity  distribution  using  the  equations  already 
developed. 
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